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PROGRESS  IN  FAVRE-REYNOLDS  STRESS  CLOSURES  FOR  COMPRESSIBLE  FLOWS 


V.  ADUMITROAIE*,  J.R.  RISTORCELLit ,  AND  D,B.  TAULBEE* 

Abstract.  A  closure  for  the  compressible  portion  of  the  pressure-strain  covariance  is  developed.  It 
is  shown  that,  within  the  context  of  a  pressure-strain  closure  assumption  linear  in  the  Reynolds  stresses, 
an  expression  for  the  pressure-dilatation  can  be  used  to  construct  a  representation  for  the  pressure-strain. 
Additional  closures  for  the  unclosed  terms  in  the  Favre-Reynolds  stress  equations  involving  the  mean  ac¬ 
celeration  are  also  constructed.  The  closures  accommodate  compressibility  corrections  depending  on  the 
magnitude  of  the  turbulent  Mach  number,  the  mean  density  gradient,  the  mean  pressure  gradient,  the  mean 
dilatation,  and,  of  course,  the  mean  velocity  gradients.  The  effects  of  the  compressibility  corrections  are 
consistent  with  current  DNS  results.  Using  the  compressible  pressure-strain  and  mean  acceleration  closures 
in  the  Favre-Reynolds  stress  equations  an  algebraic  closure  for  the  Favre-Reynolds  stresses  is  constructed. 
Noteworthy  is  the  fact  that,  in  the  absence  of  mean  velocity  gradients,  the  mean  density  gradient  produces 
Favre-Reynolds  stresses  in  accelerating  mean  flows.  Computations  of  the  mixing  layer  using  the  compress¬ 
ible  closures  developed  are  described.  Full  Reynolds  stress  closure  and  two-equation  algebraic  models  are 
compared  to  laboratory  data.  The  mixing  layer  configuration  computations  are  compared  to  laboratory 
data;  since  the  laboratory  data  for  the  turbulence  stresses  is  inconsistent,  this  comparison  is  inconclusive. 
Comparisons  for  the  spread  rate  reduction  indicate  a  sizable  decrease  in  the  mixing  layer  growth  rate. 

Keywords:  compressible  turbulence,  compressible  pressure-strain,  algebraic  Reynolds  stress  models 

Subject  classification:  Fluid  Mechanics 

1.  Introduction.  There  has  been  a  resurgence  of  interest  in  effects  of  compressibility  on  turbulent 
flows  related  to  the  design  of  high-speed/high- altitude  engines.  Although  experimental  and  numerical  infor¬ 
mation  is  growing  (for  reviews  see  [1,  2,  3])  rational  theoretical  and  modeling  efforts  are  in  a  preliminary 
stage  of  development.  This  is  consistent  with  the  fact  that  what  is  understood  of  the  relative  importance  of 
many  of  the  different  physical  effects  of  compressibility  is  very  much  in  flux  —  changing  as  more  numerical 
simulation  information  becomes  available.  While  several  issues  regarding  the  effects  of  the  compressibility  of 
the  turbulent  fluctuations  have  been  recognized,  progress  in  incorporating  the  physics  responsible  for  these 
new  compressibility  effects  in  single-point  turbulence  closures  has  been  slow.  One  of  the  major  impediments 
to  progress  has  been  the  absence  of  a  procedure  that  would  allow  for  the  inclusion  of  the  effects  of  com¬ 
pressibility  on  the  pressure-strain  covariance  appearing  in  the  second-order  moment  equations.  A  method  of 
including  compressibility  effects  as  they  appear  in  the  pressure-strain  covariance  and  also  the  variable  inertia 
effects  are  the  subject  of  the  present  study. 
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Several  earlier  studies  have  obtained  closures  for  diverse  effects  of  compressibility.  Researchers  have,  in 
general,  exploited  a  decomposition  of  the  compressible  field  into  solenoidal  and  dilatational  parts.  This  has 
been  done  using  a  dimensional  analysis  in  physical  space  [4]  or  in  Fourier  space  [5],  asymptotic  analysis  [6], 
rapid-distortion  theory  [7],  and  a  singular  perturbation  method  [8].  All  these  approaches  have  generated 
models  for  new  compressible  scalar  terms,  the  pressure-dilatation  and  the  dilatational  dissipation,  that 
appear  in  the  kinetic  energy  equation  for  high-speed  flows.  Such  an  approach  to  compressible  turbulence 
modeling  has  been  called,  very  sensibly,  an  “energetic”  approach  to  the  effects  of  compressibility,  Simone 
et  al  [9].  The  models  resulting  from  the  “energetic”  approach  to  compressibility  have  been  applied  as 
compressibility  corrections  to  the  standard  fc  -  e  model,  [10],  and  their  generalizations,  [5],  or  to  standard 
incompressible  second-order  moment  [Reynolds  stress]  closures  [11,  12,  13].  Such  a  procedure  is,  of  comse, 
an  implicit  statement  that  compressible  effects  do  not  manifest  themselves  in  either  the  pressure-strain  or 
in  the  dissipation  of  enstrophy;  the  effects  of  compressibility  occur  only  in  those  terms  explicitly  linked  to 
the  new  terms  involving  the  dilatational  field.  Thus,  energetic  approaches  to  the  problem  of  compressible 
turbulence,  as  pointed  out  in  [9],  are  incomplete. 

At  one  time,  the  pressure- dilatation  and  the  compressible  dissipation,  on  which  modeling  effort  has  been 
expended,  were  believed  to  be  the  primary  physical  effects  contributing  to  the  reduced  growth  rate  of  the 
compressible  mixing  layer.  Recent  studies  [14,  15]  have  demonstrated  that  the  dilatational  effects  on  the 
mixing  layer  are,  in  fact,  much  smaller  than  once  believed.  In  addition,  more  recent  direct  numerical 
simulations  (DNS)  suggest  that  the  pressure-dilatation  covariance  is  nominally  more  important  than  the 
compressible  dissipation,  contrary  to  early  proposals,  [4,  6].  The  pressure- dilatation  does  not,  however, 
account  for  the  reduced  growth  of  the  mixing  layer.  It  appears,  as  suggested  in  [16],  that  the  phenomena 
responsible  for  the  reduced  growth  rate  of  the  turbulent  shear  flows  is  due  to  the  reduction  in  the  Reynolds 
shear  stress  anisotropy;  this  effect  is  thought  to  be  due  to  the  effects  of  compressibility  on  the  pressure-strain 
covariance.  This  viewpoint  is  consistent  with  the  earlier  numerical  studies  of  [17].  The  articles  [13]  and  the 
later  [18]  studied  Reynolds  stress  closures  in  the  context  of  the  DNS  of  the  homogeneous  shear.  The  study 
showed  that  the  inclusion  of  the  [then]  current  compressible  dissipation  and  pressure-dilatation  models  in  the 
Reynolds  stress  turbulence  closures  led  to  improved  predictions  for  the  turbulent  kinetic  energy.  However, 
there  were  no  changes  in  the  anisotropy  consistent  with  that  seen  in  DNS.  The  authors  concluded  that 
the  [then]  current  models  were  deficient  primarily  in  the  modeling  of  the  pressure-strain  covariance  which 
controls  the  level  of  Reynolds  stress  anisotropy. 

It  should  be  noted  that  the  improved  agreement  for  the  time  evolution  of  the  kinetic  energy,  [13,  18],  when 
using  the  compressible  models  of  that  time  cannot  be  taken  to  indicate  that  such  flows  were  rationally 
predicted.  The  models  current  at  that  time,  designed  with  erroneous  assumptions  regarding  the  importance 
of  the  dilatational  effects,  were  providing  the  right  amount  of  dissipation  but  by  mechanisms  that  did  not 
reflect  actual  flow  physics.  This  has  since  been  substantiated  numerically  in  the  studies  [9,  14,  15,  16];  all  of 
which  indicate  the  lack  of  significance  of  both  the  pressure-dilatation  and  the  compressible  dissipation.  The 
fact  that  the  compressible  dissipation  and  the  pressure-dilatation  are  nominal  effects  is  also  consistent  with 
the  analytical  development  of  [8] .  In  [8]  the  pressure-dilatation  is  shown  to  vanish  as  turbulence  approaches 
equilibrium;  the  simulations  mentioned  are  quasi-equilibrium  flows.  It  should  be  mentioned  that  in  the 
homogeneous  shear  simulations,  arguably  the  most  non-equilibrium  of  the  benchmark  flows,  the  pressure- 
dilatation  is  small  but  non-negligible;  it  is  some  5-10%  of  the  dissipation.  Our  position  is  the  same  as 
the  position  of  Vreman  et  al.  [14]  as  pithily  summarized  in  their  conclusion;  to  paraphrase,  turbulence 
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models  constructed  using  the  dilatational  dissipation  or  pressure-dilatation  to  explain  the  suppression  of  the 
turbulence  [in  the  mixing  layer]  do  not  appear  to  be  reflecting  the  correct  physics. 

Thus,  if  it  is  assumed  that  the  primary  source  of  the  reduced  mixing  rate  in  the  mixing  layer  is  due  to  the 
reduction  in  the  shear  stress  anisotropy  as  indicated  by  DNS  then  a  closure  of  the  second-moment  or  Reynolds 
stress  equations  is  required:  the  pressure-strain  covariance  appearing  in  the  second-moment  equations  is  the 
only  possible  mechanism  for  such  behavior.  This  is  a  substantially  more  difficult  problem  than  that  treated 
using  the  energetic  approach;  the  quantities  requiring  closure  are  no  longer  scalars  but  second-order  tensors. 
As  one  might  expect,  true  compressible  second-order  modeling  attempts  are  few,  [7,  19]. 

This  article  describes  the  development  of  a  closure  for  the  compressible  aspects  of  pressure-strain  covariance 
appearing  in  the  Favre-Reynolds  stress  equations.  [As  is  clear  from  the  material,  the  phrase  Reynolds 
stresses  will  be  used  to  refer  to  the  Favre-Reynolds  stresses.]  Closures  for  the  unclosed  terms  involving  the 
mean  acceleration  are  also  constructed.  In  as  much  as  most  engineering  calculations  [in  this  country]  are 
done  with  lower  order  closures  this  article  therefore  also  includes  the  additional  development  of  the  second- 
order  moment  closure  into  an  algebraic  Reynolds  stress  closure,  suitable  for  ffows  in  structural  equilibrium, 
following  established  procedures  [20,  21,  22,  23,  24,  25,  26,  27,  28]. 

The  algebraic  Reynolds  stress  expression  is  noteworthy  for  the  fact  that  it  indicates  that,  even  in  the  absence 
of  mean  deformation,  the  mean  density  gradient  is  a  source  of  turbulence  stresses  in  accelerating  mean  ffows. 
As  a  consequence,  for  ffows  with  large  arbitrary  mean  density  and  pressure  gradients  an  eddy  viscosity 
representation  for  the  Reynolds  stresses  is,  from  first  principles,  inappropriate.  Qualifications  regarding  this 
statement  are  further  discussed  in  §4. 

In  the  next  section,  §2,  the  Favre  averaged  nondimensional  form  of  the  governing  equations  are  given.  Both 
first  and  second-order  moment  equations  for  a  compressible  medium,  with  no  combustion,  are  given.  Issues 
related  to  moment  closures  for  compressible  turbulence  are  also  outlined. 

The  development  of  closures  for  the  effects  of  compressibility  in  the  Reynolds  stress  equations  is  described 
in  §3.  It  is  shown  that,  within  the  context  of  a  pressure-strain  closure  linear  in  the  Reynolds  stresses,  an 
expression  for  the  pressure-dilatation  covariance  can  be  used  to  construct  the  off-diagonal  components  of 
the  pressure-strain  covariance.  In  this  way  the  results  of  previous  so-called  energetic  approaches,  [9],  to  the 
effects  of  compressibility  can  be  built  into  the  deviatoric  portion  of  an  expression  for  the  pressure-strain. 
It  is  hoped  that  such  a  procedure  would  allow  one  to  avoid  the  development  of  a  whole  new  theory  and 
methodology  for  a  compressible  pressure-strain  closure. 

Closures  for  the  effects  of  the  mean  acceleration,  which  involve  the  mass  flux,  are  also  developed.  The  new 
closures  account  for  the  influence  of  the  turbulent  Mach  number,  and  the  mean  density  and  pressure  gradients 
through  a  new  quantity,  the  baroclinic  dyad.  The  effects  of  the  bulk  dilatation  are  also  included.  Section 
3  is  concluded  with  a  summary  of  models  for  the  compressible  dissipation;  given  the  acknowledged  lack 
of  importance  of  the  compressible  dissipation  in  weakly  compressible  aerodynamic  tiubulence,  as  discussed 
above  and  in  §3,  no  development  of  models  for  the  compressible  dissipation  is  pursued. 

Starting  with  the  closure  for  the  second-order  moment  equations  developed  in  §3,  §4  develops  an  algebraic  clo¬ 
sure  for  the  Reynolds  stresses.  The  physics  of  compressibility  as  captured  by  the  full  second-moment  closure 
are  therefore  built  into  the  simpler  and  more  widely  used  two-equation  k  -  e  platform.  The  tensor  polyno¬ 
mial  representation  techniques  employed  produce  both  two-dimensional  and  three-dimensional  versions  of  an 
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algebraic  turbulence  stress  closure.  Given  the  complexity  of  the  three-dimensional  algebraic  closure  and  the 
current  status  of  single-point  turbulence  models  for  three-dimensional  flow  only  the  two-dimensional  model 
is  developed  into  a  working  closure. 

Section  5  focuses  on  a  numerical  investigation  of  the  closures.  The  numerical  method  used  to  simulate  the 
free-shear  flows  of  interest  is  sketched.  The  theory  and  results  presented  in  earlier  sections  are  implemented 
in  simulations  in  the  mixing  layer. 

Simulations  using  second-order  moment  (SOM)  closures  as  well  as  the  algebraic  models  are  conducted.  The 
numerical  experiments  are  constructed  with  the  intention  of  investigating  several  very  different  issues  of 
relevance  to  the  prediction  of  compressible  turbulent  flows  for  engineering  purposes.  Foremost  in  importance 
is  the  ability  of  an  algebraic  stress  model  -  that  includes  the  effects  of  compressibility  -  to  reflect  the  physics 
in  the  full  SOM  closures. 

The  pressure-strain  methodology  developed  for  the  SOM  equations  in  §3  are  general  and  depend  on  a  choice  of 
closures  for  the  pressure-dilatation.  As  consequence  it  follows  that  it  is  necessary  to  understand  sensitivity 
of  the  formulation  to  different  models  for  the  pressure-dilatation.  In  this  context  two  pressure-dilatation 
models  are  investigated. 

In  addition  to  assessing  the  sensitivity  of  the  pressure-strain  model  to  different  pressure-dilatation  models 
and  the  suitability  of  the  algebraic  stress  closure,  we  compare  the  computational  results  to  what  is  seen 
in  numerical  and  laboratory  experiments.  Of  particular  interest  is  the  well  known  effect  of  compressibility 
on  reduction  of  the  spread  rate  of  the  mixing  layer.  In  as  much  as  the  reduction  in  the  spread  rate  is 
due  to  changes  in  the  the  anisotropy  of  the  Reynolds  stresses,  as  a  function  of  compressibility,  are  also 
investigated.  Given  that  the  anisotropy  of  the  Reynolds  stresses  is  dependent  on  the  pressure-strain,  the 
effects  of  compressibility  on  the  different  components  of  the  pressure-strain  are  also  investigated. 

2.  Governing  equations.  The  problem  formulation  is  now  described.  This  includes  a  statement  of 
the  governing  equations  -  the  first  and  second- moment  equations.  Indications  of  the  modeling  issues  to  be 
addressed  in  subsequent  sections  are  also  given.  The  Favre  averaging  procedure  is  first  described. 

The  conservation  equations  for  mass,  momentum  and  energy  in  a  Favre  setting  are  now  derived.  The 
dependent  variables  used  are  the  density  p,  the  velocity  Ui  and  the  total  energy  et  =  h  —  p/p  -f  UiUi/2. 
The  fluid  is  assumed  to  be  Newtonian  fluid  satisfying  the  Stokes  relation  with  zero  bulk  viscosity  and  a 
constant  molecular  Pfandtl  number.  Although  real  gas  effects  are  of  interest  for  industrial  applications 
there  are  a  sufficient  number  of  unresolved  and  more  important  issues  associated  with  compressibility  that 
justify  limiting  the  study  to  ideal  gases.  Consequently,  the  pressure  p  is  obtained  from  the  equation  of  state 
p  =  pRT. 

The  Favre  averaging  procedure  is  now  described:  denote  by  an  over-bar  the  ensemble  (or  time)  average  and 
by  the  brackets  the  density  weighted  ensemble  (or  time)  average: 

(1)  W  =  T 

The  ensemble  average  obeys  the  following  decomposition  rules: 

(2)  X  =  = 


4 


The  Favre  average  obeys  the  following  decomposition  rules: 

(3)  X  =  {X)^X'\  (X")=0, 

The  application  of  the  above  averaging  procedure  on  the  instantaneous  transport  equations,  decomposing  the 
variables  Favre  mean  and  fluctuating  components,  produces  the  Favre  averaged  equations.  The  equations  of 
motion  are  first  rewritten  in  nondimensional  form  (with  respect  to  reference  values  taken  in  the  high  speed 
stream:  Poo,  ^oo,  ^oo  and  the  inlet  value  of  the  vorticity  thickness  5^^),  Using  these  reference  quantities, 
we  define  the  relevant  nondimensional  parameters:  the  Reynolds  number  Re  =  Pc»'^oo^w/Moo?  the  Prandtl 
number  Pr  =  Cpiioo/\  the  mean  flow  Mach  number  M  —  Uoo/VT-RToo,  7  is  Cp/Cy.  For  this  study,  Pr  =  1. 


2.1.  First  moment  equations.  The  mass  conservation  equation  reads: 


(4) 


dp  dpjuj)  _ 

dt  dxj 


and  the  conservation  of  momentum  is: 

dpjui)  dp{ui){uj)  ^  drp{u'!u'!)  _  dajiju) 

^  dt  dxj  dxj  dxi  dxj 

Note  that  the  stress  tensor  notation  is  changed  to  aji{u)  =  ^[Sij{u)  —  ^Spp{u)Sij]  =  2/jLSfj{u)/Re.  For  the 
strain  rate  Sij{u)  =  and  for  all  the  other  linear  differential  operators,  this  new  notation  is  more 

suitable  when  investigating  compressible  flows.  In  these  instances  the  two  types  of  averages,  having  different 
properties  with  respect  to  the  linear  differential  operators,  are  naturally  encountered.  Hereafter  any  tensor 
with  a  star  superscript  will  indicate  the  deviator  -  the  traceless  portion  of  that  tensor. 


Two  supplementary  hypotheses  pertaining  to  the  molecular  transport  of  momentum  are  set  forth:  the 
viscosity  fluctuations  are  unimportant  and  the  mean  viscosity  (p)  is  described  by  the  Maxwell-Rayleigh  law, 
i.e.,  it  varies  with  the  mean  temperature  as  (p)//ire/  =  =  0-76.  The  fluctuations  of  the 

viscosity  and  their  correlation  with  other  variables  are  neglected.  Within  the  current  notation  the  averaged 
stress  aji{u)  is  equal  to  (Jji{{u))  + 


The  Favre  mean  of  the  total  energy  of  the  fluid  et  =  T/[7(7  —  1)M^]  +  UjUjl2  obeys: 


(6) 


dp{et)  dp{et){uj)  _  dq^jT) 

dt  dxj  dxj 


dfpiu'fe'l)  d  , _  ^  ^ 


As  is  the  situation  with  the  momentum  equation,  the  molecular  flux  term  in  the  energy  equation  introduces 
the  other  type  of  mean  quantity  present  in  the  Favre  formulation  -  the  plain  ensemble  mean.  Using  the  new 
notation,  the  heat  flux 

^  P  dT 

’  {j  -  l)RePrM^  dxj 

has  the  resulting  averaged  expression: 

q^{T)  =  qj{{T))  +  qjiT^). 

The  remaining  terms  on  the  right-hand  side  of  the  energy  equation  can  be  expanded  by  using  the  equation 
of  state  in  the  average  sense 
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Then  the  pressure  work  term  is  of  the  form 


the  viscous  work  term  reads 


puj^p{uj)  +  -^p{u'^T"), 


Uiaji{u)  =  {ui)aji{{u))  -f  {ui)aji{u'^)  +  u'^(Xji{{u))  -f  u'\aji(u^^) 


and  the  turbulent  total  energy  flux  is 


1 

7(7  -  1)M2 


The  fluctuations  of  transported  quantities  in  turbulent  flows  are  sustained  via  interaction  between  the  mean 
gradients  and  the  turbulence.  The  fluctuations  will  be  characterized  by  length  and  time  scales  of  the  order  of 
those  of  the  mean  flow  and  which,  provided  the  Reynolds  number  is  large,  will  be  many  orders  of  magnitude 
larger  than  the  fine  scales  at  which  the  molecular  diffusion  is  important.  Thus  the  form  of  the  moment 
equations  carried  in  the  closure  development  for  the  computation  of  free  shear  flows  will  not  include  viscous 
transport  effects. 


2.2.  Second-order  moment  equations.  The  average  of  the  first  moment  of  the  Favre  decomposed 
Navier-Stokes  equations  produces  the  following  second  moment  equations: 

dp{u'lu'j)  dp{u'^u'!){uk)  _ ^ 

dt  dxk  dxk 


p{u"u'-ul)  +p'u'!5jk  -  u'!akj{u")- 


d{uj) 


(8) 


u>!aki{un\ -p{«r-^-p{<u';y-^ 


\dxj 


dxi )  ‘  dxj  ^  dxi 


■^d(Jki{{u))  ,  —dakj{{u)) 


dxk 


dxk 


It  is  necessary  to  provide  a  closure  for  two  compressible  quantities:  the  pressure-strain  correlation  and  the 
mass  flux/pressure  gradient  (the  last  three  terms  in  the  second  line).  The  molecular  diffusion  terms  are 
generally  small  in  high  Reynolds  number  flows  and  will  be  negligible  for  problems  addressed  here.  It  should 
be  pointed  out  that,  depending  on  how  the  derivation  is  done,  the  mass  flux  terms  multiply  the  acceleration 
of  the  mean  flow.  For  this  reason  the  mass  flux  terms  are  often  called  acceleration  terms:  they  appear  to  be 
important  in  accelerating  mean  flows. 


Various  tensors  are  divided  into  their  traces  and  their  deviators.  The  production  of  the  Reynolds  stresses 
becomes. 


(9) 


P*  —  P.  . 


dxk 


dXk 


djui) 

dxk 


where  P  is  the  production  of  (fc) .  The  pressure-strain  covariance  is  written  as 


(10) 


TT*_TT  2 

ny  -  ny  gP  d6ij  -p  j  gP  Sij. 


The  viscous  and  pressure  acceleration  terms  are  written,  respectively,  as 


2^^  —d(Jkj{{u)) 

dXk  ^  dXk 


2^daki{{u)) 

3  '  dxk 


(11) 


and 


(12) 


Mtj  =  Mij  - 


+7i"  — 

*  dxj  ^  dx, 


3  '^dxk 


Note  that  the  mean  pressure  gradient  appearing  in  Mij  and  the  mean  viscous  stress  appearing  in  Vij  can  be 
replaced  using  the  mean  momentum  equations.  In  which  case  they  are  written  in  terms  of  the  mean  flow 
Lagrangian  acceleration.  Both  of  these  terms  involve  the  mass  flux;  any  closure  for  these  two  acceleration 
terms  requires  a  closure  for  the  mass  flux.  The  dissipation  is  rewritten  as 


(9  \  dv/-  2 

In  turbulent  free  shear  flows  the  viscous  diffusion  part  is  overlooked  owing  to  high  Reynolds  numbers  which 
are  characteristic  to  these  flows.  The  present  analysis  can  accommodate  the  discarded  viscous  terms  when 
necessary;  the  tensor  Vy  will  be  carried  for  generality.  Furthermore,  the  turbulent  transport  terms  are 
considered  together: 

(14)  Tijk  =  -  +^Sik  +FWSjk  -  u'!cTkj{u")  -  u'!aki(u")  . 


The  Reynolds  stress  or  second-order  moment  equations  are  then  written 


(15) 


dp{u'lu'!)  dp{u>fu'J){uk)  _  ^  y  ,  P*  I  n*  I  M*  +V*  -  re* 

+l[P  +  ¥d  +  M  +  V-pi\5ij 
o 


The  equation  for  the  anisotropy.  The  single  point  anisotropy  tensor,  a,  is  defined 
(16)  a,ij  =  {u^  Uj )/ (fc)  —  ’ 

The  anisotropy  has  zero  trace.  Note  that  the  Reynolds  stress  has  been  normalized  by  (fc)  and  not  2(fc); 
in  which  case  this  version  of  the  anisotropy  is  related  to  an  older  definition  by  aij  =  2bij.  The  turbulence 
time  scale  is  r  =  (fc)  /e^  where  Is  is  the  solenoidal  dissipation.  The  second  invariants  of  the  mean  strain  and 
rotation  are  denoted  by  =  Sfj{{u))Sjj^{{u))  and  o;^  =  ^ij{{n))Q>ij{{u)). 

For  subsequent  manipulations  used  to  derive  an  algebraic  stress  closure  in  §4  we  replace  the  second-order 
transport  equations  with  an  equation  for  (fc)  and  the  anisotropy  tensor.  The  following  form  of  the  transport 
equation  for  the  anisotropy  tensor,  [21],  will  be  used: 


(17) 


rap- 


Dt 


1 

dTijk 

{uy;)  dTk] 

Oij 

\dTk 

(fe) 

dxk 

(k)  dxk_ 

(k) 

dxk 

dxk_ 

— 

ij. 


4  -  2  +  (2  - 

p  Cs  cr  Dt 


P  Cs 


where  D/Dt  indicates  the  mean  Lagrangian  derivative,  DfDt=^  ^  -h 
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The  kinetic  energy  equation.  The  equation  for  the  turbulent  kinetic  energy  (fc)  =  {u'-u'-)/2  is 

d  r_,  .  -r-^  — - _/  //  ,nd(ui) 


dp{k)  dp{k){uj)  __ 
dt  dxj 


(18) 


—dc^jiiju)) 


.du'J 


du'' 


A-  n'rri.  _  TTi  I  _  ff..(uff)irL 


The  additional  terms  reflecting  the  compressible  nature  of  the  turbulence  are  on  the  last  line;  they  are, 
respectively,  the  pressure-dilatation  covariance,  the  mass  flux  -  mean  acceleration,  and  the  compressible 
dissipation.  Note  that  neither  the  pressure-strain  nor  the  mass  flux/ pressure  acceleration  terms,  Aiij ,  appear 
in  the  (fc)  equation.  The  effects  of  the  pressure-strain  and  the  Mij,  appears  only  in  the  Reynolds  stress 
equations.  Any  classical  two-equation  turbulence  model  cannot,  as  a  consequence,  account  for  the  physics 
associated  with  these  two  unknown  covariances  that  lead  to  changes  in  the  Reynolds  stress  structure.  Early 
approaches  accounting  for  the  effects  of  compressibility  have  focussed  only  on  the  effects  of  compressibility 
as  they  occur  in  the  energy  equation.  This  so-called  “energetic”  approach,  to  use  the  phrase  of  Simone  et  al. 
[9],  for  the  effects  of  compressibility  misses  the  changes  in  the  Reynolds  stresses  due  to  compressibility  and 
inertia.  It  is  for  this  reason  that,  in  §3,  the  second-order  equations  are  closed.  An  algebraic  Reynolds  stress 
closure  is  then  derived  in  §4.  In  this  way  the  structural  effects  of  compressibility,  for  a  certain  class  of  flows, 
can  be  accounted  for  in  the  context  of  a  two-equation  single-point  closure. 


The  dissipation  equation.  The  dissipation  appearing  in  the  kinetic  energy  equation  is  typically  written, 
for  locally  homogeneous  flows. 


(19) 


.du>/ 


=  P€  =  p{es  +  ec)  =  2pQ,ij{u")Q,ij{u")  -h 


[4,  6,  29],  with  fly(u)  =  -  §^)  as  the  rotation  rate.  It  is  customary  to  compute  the  solenoidal 

dissipation  from  the  incompressible  k  —  e  model  extended  to  variable  density  flows: 


(20) 


dp  ,  dp  es{uj)  ^ ^  T  .  „  _  p  dls 

dt  dxj  dxj  _  ^  ^  Re  dxj 


djuj) 

dXj 


where  Cei  =  1.44  and  Ceg  =  1-92. 


3.  Compressible  closures  for  the  second-order  moment  equations.  This  section  describes  the 
development  of  closures  for  the  unknown  terms  in  the  second-order  equations.  The  pressure-strain  covari¬ 
ance  appearing  in  the  Reynolds  stress  equations  is  closed  using  a  linear  tensor  polynomial  representation 
in  the  Reynolds  stress  following  well  established  procedures  similar  to  Launder  et  al  [30].  In  recogni¬ 
tion  of  the  minor  role  viscosity  plays  in  high  Reynolds  number  free  shear  flows,  only  pressure-strain  and 
pressure-acceleration  closures  are  addressed.  The  pressure- acceleration  terms  are  closed  using  a  leading  or¬ 
der  [isotropic]  gradient  transport  expression  for  the  mass  flux.  The  following  section,  §4,  uses  the  results  of 
the  present  section  to  obtain  an  algebraic  closure  for  the  Reynolds  stresses. 

3.1.  A  closure  for  the  pressure-strain  covariance.  In  incompressible  turbulence  the  closure  of  the 
pressure-strain  covariance  by  a  tensor  polynomial  linear  in  the  Reynolds  stresses  is  used  with  notable  success 
for  simple  two-dimensional  mean  flows.  The  procedure  is  standard;  a  popular  early  reference  is  Launder 
et  al[S0].  An  updated  version  of  the  linear  [30]  rapid  pressure-strain  modeling  is  given  in  Speziale  et  al 
[31].  As  the  ultimate  goal  is  to  devise  an  algebraic  Reynolds  stress  model  only  pressure-strain  models  linear 
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in  the  Reynolds  stresses  will  be  considered.  A  discussion  of  nonlinear  rapid-pressure-strain  model  and  the 
need  for  realizability  can  be  found  in  [32].  Additional  discussion  of  the  physical  assumptions  underlying  such 
methodologies  can  be  found  in  [33],  [34]  and  [35]. 

The  compressible  correction  to  the  pressure-strain  covariance  representation  is  obtained  using  established 
linear  procedures.  As  in  the  incompressible  situation,  [30,  34,  32],  the  pressure  strain-covariance  closure  can 
be  written  as: 


(21)  Yiij  —  P  —  A.ij  -h  2p  ^piqj  +  [5pg((u))  -h  fipg((w))]* 

Such  an  expression  is  possible  if  the  supplementary  compressible  terms  appearing  in  the  Poisson  equation 
for  the  pressure  (see  Appendix)  are,  in  the  weakly  compressible  limit,  of  higher  order.  This  will  be  the  case 
if  the  evanescent  wave  portion  of  the  initial  value  problem  is  not  important  as  is  the  case  for  aerodynamic 
applications,  [36].  The  tensors  in  the  above  decomposition  are  modeled  as: 


(22) 


(23) 


—  C7ip  cdij  “j“  '>/App  ^ 

=  Oi\SqiSpj  -|-  Oi2(^Spq^ij  “f"  ^qj^pi^  4“  01^pj^qi 
“h  /^2(^pqr®ij  d”  ^pi^qj  "h  ^ij^pq  d~  ^jq^pi^  "h  f^S^qi^pj 


Here  Ipjqi,  following  precedent,  is  linear  in  the  Reynolds  stresses  satisfying  necessary  symmetry  requirements. 
There  are  five  unknowns.  To  determine  the  coefficients  in  Ipjqi  additional  constraints  are  required.  As  is 
the  usual  procedure,  the  normalization  constraint,  [30,  32],  requires 

(24)  ^ppqi  —  [(3ai  d-  2a2)Sqi  +  (3/3i  +  4^2)<^gi](fc)  —  i'^q'^i) 


which  can  be  satisfied  if 


3q:i  d"  2qj2  —  2/3 

3/3i+ 4/32  =  1. 

There  are  now  three  unknowns  and  additional  information  is  required  to  obtain  them. 

In  incompressible  turbulence  modeling  the  trace  of  the  pressure-strain  is  zero,  Ipkqk  —  0,  and  this  provides 
the  additional  information  to  determine  the  unknowns.  In  compressible  turbulence  the  trace  of  the  pressure- 
strain  is  the  pressure-dilatation  and  the  so-called  continuity  constraint  becomes 

(25)  *A.kk  d"  4/?  Ipkqk  [^pq  d*  ^pq]  —  2p'd. 

As  the  right  hand  side,  pd,  is  known  from  earlier  energetic  approaches  to  the  compressible  turbulence  model¬ 
ing  problem,  the  continuity  constraint,  equation  (25),  becomes  a  constraint  that  determines  the  coefficients 
in  the  pressure-strain  closure.  Note  that  as  pd  vanishes  with  turbulent  Mach  number  the  incompressible 
limit,  Akk  +  Ipkqk[Spq  +  fJpg]  =  0,  is  recovered  for  vanishing  compressibility. 

As  a  consequence  of  the  continuity  constraint  a  certain  combination  of  coefficients  appears  in  the  final  model. 
These  are  readily  defined  by  a  portion  of  (25)  such  that 

(26)  Ipiqi  f  (^)  —  (^1  d"  4ck2)^pqf  d”  (/3l  d"  6^2  d~  fd^CLpq  —  d\Spq  -(-  d2Clpq, 
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The  final  model  for  the  compressible  portions  of  the  pressure-strain  is  then  written  in  terms  of  di  and 
^2-  The  values  of  di  and  ^2,  as  will  be  indicated  shortly,  are  then  determined  by  the  expression  for  the 
pressure-dilatation.  The  values  of  the  at  can  be  related  to  the  df. 

ai  = -di/5  +  4/15,  a2  =  3di/10-l/15 
/?i  =  (15  +  6C2)/33,  /?2  =  -(2  +  3C2)/22 
p3  =  d2 


and  Ca  =  (5  --  QCs)/!!,  C4  -  (1  +  TCs)/!!. 

Application  of  the  normalization  and  continuity  constraints  then  allows  the  linear  pressure-strain  model  to 
be  written  as 


(27) 


n*-  - 


[1  —  C3  +  2^2] 


p  =  -Cip  eaij  +  p{k) 


+ 1*)  S5.((«))+ 

2 


[1  —  (^4  2^2]  ((Ur))  ^^ip({w))^pj]  “1“  0  d2<S^pp((‘U))(2ij 


which  follows  from  (21)  after  subtracting  the  pressure-dilatation  from  the  left  side  and  the  pressure-dilatation 
model  from  the  right  side  of  (21).  All  terms  involving  the  di  represent  corrections  due  to  the  compressibility 
of  the  fluctuations;  the  di  vanish  as  the  turbulent  Mach  number  vanishes.  The  terms  involving  the  Ci 
come  from  the  incompressible  pressure-strain  model.  In  which  case  the  choice  of  the  Ci  allow  one’s  favorite 
incompressible  pressure-strain  model  to  be  used.  The  above  expression  for  the  pressure-strain  is  used  in 
calculations  presented  in  this  paper.  The  calculations  in  this  article  are  of  standard  test  flows  and  are  very 
simple  flows.  In  more  complex  flows  realizable  issues  may  need  to  be  incorporated.  A  realizable  form  of  the 
model  is  given  in  an  Appendix;  more  details  can  be  found  in  Adumitroaie  [37]. 


3.1.1.  Commentary.  It  has  been  seen  that  the  knowledge  of  one  invariant  of  the  pressure-strain, 
the  pressure-dilatation,  and  the  assumption  of  a  form  linear  in  the  Reynolds  stress  allows  one  to  obtain  a 
model  for  the  deviatoric  components  of  the  pressure-strain  tensor.  In  this  way  an  independent  new  theory 
for  the  compressible  pressure-strain  can  be  avoided  by  using  the  results  of  developments  for  the  scalar 
pressure- dilatation.  The  results  of  previous  so-called  energetic  approaches,  Simone  et  al.  [9],  to  the  effects 
of  compressibility  is  built  into  the  deviatoric  portion  of  an  expression  for  the  pressure-strain. 


This  point  merits  consideration  from  another  point  of  view.  Consider,  for  the  moment,  the  following  partition 
of  the  pressure-strain: 


(28) 


psij  =  ps^j  +  -pd6ij. 


For  this  subsection,  the  primes  on  p  have  been  dropped.  Here  by  s*j  we  now  mean  the  deviatoric  portion 
of  the  fluctuating  strain,  Sjj  =  0,  which  contains  solenoidal  and  compressible  contributions  s*j  =  s'-j  -h 
where,  of  course,  .  This  is  done  for  ease  of  presentation;  in  the  nomenclature  of  §2  these  quantities 

would,  of  course,  be  represented  by  S*j{u")  a  precision  unnecessary  for  the  present  discussion.  Thus  one  can 
write 


(29) 


pSij  = 


+  P^ii 


+ 


-  pd  6ij . 


10 


The  term  is  closed  using  standard  incompressible  pressure-strain  closures.  The  pressure-dilatation 

is  closed  using  models  already  in  the  literature,  see  §3.2  below.  In  §3.1  an  expression  for  has  been 
obtained. 


As  has  already  been  discussed  in  energetic  approaches,  the  pressure-dilatation  cannot  account  for  the  suppres¬ 
sion  of  the  turbulence  by  the  reduction  in  the  shear  stress  as  is  seen  in  compressible  flows.  A  straightforward 
extension  of  the  energetic  approach  to  the  second-order  closure  level  implies 


(30) 


pSij  =  ps*i.  + 


Such  an  expression,  as  is  consistent  with  [13,  18],  is  likewise  unsuccessful  in  reducing  the  turbulence  shear 
stress.  As  will  be  discussed  further  in  §5  the  pressure-dilatation  is  a  small  quantity  and  has  a  nominal 
effect,  for  a  «  1,  on  the  turbulence  stresses  and  energy  when  it  is  included  in  the  spherical  portion  of  the 
pressure-strain.  Our  procedure,  using  the  pressure-dilatation,  produces  an  expression  for  the  compressible 
contribution  to  the  deviatoric  portions  of  the  pressure  strain,  ps*j  .  In  fact,  as  was  born  out  by  computational 
experiments,  had  the  expression 


(31) 


pSij  =  ps\ 


-I  , 

f  +Ki 


[without  pd  on  the  diagonal]  been  used  in  the  numerical  investigations  reported  in  §5,  our  results  would  not 
have  changed  much,  for  a  <  1. 


3.1.2.  Pressure-dilatation  models.  To  obtain  the  final  form  of  the  pressure-strain  an  expression  for 
^  is  required.  There  are  some  choices.  Some  proposals  for  pd  require  the  solution  of  transport  equations 
for  the  density  variance  [5]  or  pressure  variance  [7].  There  are  two  pressure-dilatation  models  [8,  38]  that  do 
not  require  separate  equations.  The  model  of  Sarkar  [38]  is; 


(32)  ln„  =  p-d  =  -  g? 

where  xi  =  0-15,  X2  =  0.2,  and  Xs  —  0  (still  to  be  determined  by  the  author).  Note  that  the  model  as  we 
have  grouped  the  terms  makes  it  appear  as  if  it  is  singular  in  Mt\  it  is  most  definitely  not  singular.  The 
model  of  Ristorcelli  [8]  is 


p'd  =  -xM‘^[P  -  p  e^-Tk-  -  1){Pt  +  p  e  +  Tt)] 


(33) 


,£»(3a2-h5a;2) 
Dt 


where 

/  _  _ ^pd _ 

1  +  2IpdMi  +  llpdMtiii  -  1)  ’  ^  ~  1  +  2IpdM?  -t-  llpdMhi'i  -  1)  ’ 

(34)  Ipd  =  g/i  +  Ipd  [3cr^  -t- buF\  ,  ^  ^3) 


Here  a  is  the  proportionality  constant  in  the  Kolmogorov  scaling;  we  return  to  this  in  more  detail  below. 
These  two  models,  as  well  as  a  third  (the  Aupoix  model),  have  been  discussed  in  [39]. 


The  di  are  determined  by  the  pressure-dilatation  model.  For  the  Sarkar  model,  one  requires  that 


(35) 


di 


SxaMf  ,  ^  XiMt 
3  ’  ^  2 
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(36) 


^  =  X2Mfp  ts. 


For  the  Ristorcelli  model  one  requires  that 


(37) 


,  xMf  ,  xMi 
d,  =  —d2  =  — 


(38)  ^  =  xMtlp  e  +  -  1){Pt  +  p£)]-  p{k)M^x'^^^^-^^^- 

Transport  terms  have  been  neglected;  this  is  necessary  to  obtain  an  algebraic  closure  [which  requires  ho¬ 
mogeneity].  The  compressibility  of  the  turbulence  will  only  be  important  where  the  turbulence  energy  is 
large  [and  thus  also  Mt],  which  is  typically  in  regions  of  large  production  and  where  transport  is  not  as 
important.  Which  is  to  say  that  transport  is  only  of  importance  in  the  peripheral  regions  of  simple  flows  - 
regions  where  the  turbulence  intensity  is  low,  production  low,  and  [therefore]  the  fluctuations  are  essentially 
incompressible. 


3.2.  The  mass  flux  and  pressure-acceleration  closures.  The  models  for  mass  fluxes  have  not 
undergone  much  development  and  their  importance  in  a  general  flow  is  not  fully  understood.  Judging  from 
the  equations  one  can  infer  that  the  mass  flux  will  be  important  in  accelerating  flows  with  large  density 
gradients.  There  has  been  some  research  on  the  mass  flux:  [5]  solved  its  transport  equations  and  Ristorcelli 
[40]  showed  how  an  algebraic  model  for  the  mass  flux  can  be  obtained  from  its  transport  equation: 


(39) 


P  u'i 


r  . 


dxk  dx 


'j  J 


where  Tu  =  Mtr/[1  (P  ^)  “  !)]♦  Here  i/q?  and  1/2  are  the  coefficients  retrieved  from  the  inversion  of 

the  matrix  Gij  =  Sij  -hTu^^:  I'o  =  Ig-\- ^2  =  il  +  lG  +  IlG  +  IIlG)~^,  the 

Roman  numbers  representing  the  invariants  of  G.  Note  that  the  leading  order  term  is  a  gradient  transport 
model.  For  an  isotropic  turbulence,  an  eddy  viscosity  formulation  is  possible.  For  the  present  set  of  simple 
benchmark  shear  flows,  one  does  not  expect  the  mass  flux  terms  to  be  very  important  -  the  mean  flow  does 
not  accelerate  much.  We  shall  for  the  sake  of  simplicity  use  the  eddy  viscosity  form  of  the  mass  flux  - 
the  lowest  order  contribution  in  the  polynomial  given  above.  This  is  consistent  with  the  gradient  transport 
model  used  in  [11]. 


The  baroclinic  dyad,  a  tensor  product  formed  from  the  mean  density  and  pressure  gradients,  is  defined  as 

dp  dp 


(40) 


dxi  dxj 


and  the  mass  flux/ acceleration  terms  in  the  second-moment  equations  can  then  be  written,  using  the  leading 
order  portion  of  39,  as 

(41)  M*J  =  -—TuMkmj  +  Tin  -  \n„6ij)  =  -i|r„i/o(fe)i?r,. 

where  =  TZij  +  IZji  -  ^TZppSij. 
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3.3.  The  compressible  dissipation.  There  are  several  models  available  for  the  compressible  dis¬ 
sipation,  [4,  5,  6,  8].  Many  of  these  models  reflect  certain  assumptions  regarding  the  importance  of  the 
compressible  dissipation  observed  in  early  DNS  of  compressible  turbulence.  The  compressible  dissipation 
has  since  been  found  to  be  less  important  than  originally  believed.  In  fact,  the  low  turbulent  Mach  number 
asymptotics  of  [8]  indicate  it  varies  inversely  with  the  Reynolds  number  and  as  a  consequence  is  negligible  in 
engineering  flows,  though  not  in  low  Reynolds  number  DNS.  Compressible  dissipation  models  are  nonetheless 
included  for  completeness. 

As  the  Taulbee  and  Van  Osdol  [5]  compressible  dissipation  model  requires  additional  transport  equations, 
in  the  spirit  of  computational  simplicity,  that  model  will  not  be  used.  The  closure  proposed  by  Sarkar  et  al 
[6]  is  based  on  ideas  from  linear  acoustics  and  appears  related  to  the  initial  value  problem,  [36].  It  can  be 
written  as 

(42)  ec  =  3asesMf 

with  tts  =  1.0  from  DNS  of  decaying  compressible  turbulence;  unfortunately  this  arrangement  was  deemed 
to  lack  universality  [19].  Zeman  [4]  provides  a  model  on  the  grounds  that  eddy  shocklets  occur  in  high  speed 
flow  and  relating  this  assumption  to  the  dilatational  dissipation: 

(43)  Ce  =  0.75(1  -  exp{-[(^(l  +7)Mt  -  0.1)/0.6]2})e, 

It  has  been  shown  by  Blaisdell  et  al.  [41]  that  Zeman’s  model  gives  incorrect  scaling  between  Cs  and  Cq. 
Besides,  the  exponential  dependence  on  Mf  delivers  a  steeper  growth  rate  reduction  compared  to  other 
models. 

As  has  been  mentioned  in  more  recent  DNS  studies,  [9,  14,  15,  42],  have  demonstrated  that  dilatational 
covariance  closures  with  a  Mf  scaling  predict  effects  of  compressibility  when  they  are,  in  fact,  very  small. 
In  an  asymptotic  analysis,  Ristorcelli  [8]  has  found  that  the  compressible  dissipation  has  an  dependence 
and  is  inversely  proportional  to  the  turbulent  Reynolds  number: 


The  parameters  are  If  =  0.3,  /|  =  13.768,  /|  =  2.623,  /[  =  1.392,  =  3  and  =  0.4-4  is  the 

Kolmogorov  scaling  coefficient.  Also,  Mt  denotes  the  turbulent  Mach  number,  that  is  Mf  =  |(A;)/c^  and  Rt 
the  turbulent  Reynolds  number  Rt  =  p4{k)‘^/{9eati)Re.  The  local  speed  of  sound  is  given  by  =  T/M^. 
For  high  Reynolds  number  flows,  in  the  absence  of  wall  effects,  viscous  diffusion  is  negligible. 

At  this  point  the  Reynolds  stress  equations  have  been  closed  and  it  is  possible  to  compute  the  flow  using  a 
second-order  closure.  This  is  the  subject  of  §5.  However,  there  is  a  wide  class  of  flows  of  engineering  interest 
that  can  be  treated  using  simple  algebraic  Reynolds  stress  models.  An  algebraic  closure  for  the  Reynolds 
stresses  is  therefore  first  developed. 

4.  A  compressible  algebraic  Reynolds  stress  model.  A  quasi-explicit  algebraic  model  for  the 
Reynolds  stresses  is  now  derived.  An  algebraic  Reynolds  stress  model  comes  from  the  fixed  point  solution 
of  the  evolution  equations  for  the  anisotropy  tensor.  These  equations  can  be  thought  of  as  describing  a 
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turbulence  in  a  state  of  structural  equilibrium:  the  fixed  point  solution  corresponds  to  an  exact  solution  of 
the  Reynolds  stress  equations. 

Several  permutations  of  quasi-explicit  algebraic  Reynolds  stress  expressions  exist:[20,  21,  22,  23,  24,  25,  26, 
27,  28].  The  qualifier  “quasi-explicit”  is  used  to  indicate  that,  as  the  fixed  point  equations  are  nonlinear, 
the  solution  is  given  implicitly.  A  notable  exception  is  the  recent  explicit  algebraic  model  of  Girimaji  [28] 
who  has  found  the  exact  nonlinear  solution  to  the  fixed  point  equations.  The  inception  of  our  work  predates 
[28]  and  our  procedure  follows  precedents  set  by  [20,  21,  26].  The  polynomial  representation  methods  will 
be  used  to  obtain  two-dimensional  and  three-dimensional  versions  of  the  algebraic  turbulent  stress  models. 


We  use  the  ansatz  introduced  by  Taulbee  [21],  to  obtain  the  fixed  point  equation  corresponding  to  the 
differential  equation:  the  Lagrangian  derivative  in  (17), 


(45) 


R  pi]  =  0, 

Dt  '■ra' 


is  set  to  zero,  which  allows  relaxation  effects  to  be  built  into  the  Reynolds  stress  model.  Equation  (45)  allows 
a  relaxation  of  the  anisotropy  to  its  equilibrium  value  at  the  same  rate  that  the  relative  strain  relaxes  to 
its  equilibrium  value,  [21].  The  combination  of  transport  terms  in  (17)  is  set  to  zero  following  established 
algebraic  stress  modeling  procedures: 


(46) 


1 

dTijfc 

(«)  dTk] 

aij 

W) 

_  dxk 

(fc)  dxk_ 

(k) 

dxk 

dxk^ 

Applying  these  approximations  results  in  a  quasi-linear  tensor  expression  for  the  anisotropy. 


If  6i  =  I  -  |di,  62  =  63  =  Cs  -  2^2,  b4  =  C4  +  2d2,  and 


(47) 


ai  +  a.  -  2 + {2  -  c..)=L  + +  ^,1  -  2d,)s„m  + 

€s  p  Cs  (X  JJt  6  p  es 


-1 


the  algebraic  fixed  point  form  of  the  Reynolds  stress  anisotropy  equation  is  written 

2 


(48) 


a  =  -gr 


biS*  +  62R.*  +  bs  (aS*  +  S*a  -  ^{aS*}5  )  -  btiaU  -  fia) 


0-' 


where  the  curly  braces  signify  the  trace.  Prom  this  expression  it  is  seen  that  the  anisotropy  tensor 


(49) 


aij  — 


is  dependent  on  three  second  order  tensors,  two  symmetric  and  one  skew-symmetric. 


4.1.  A  three-dimensional  algebraic  Reynolds  stress  model.  Standard  representation  theory 
methods  can  be  applied  to  obtain  the  solution  aij  -  ft,  R*)  to  (48).  In  contrast  to  the  incompressible 

case  the  solution  of  (48)  is  now  much  more  difficult  inasmuch  as  the  procedure  now  involves  an  additional 
tensor.  Following  standard  methods  the  solution  can  be  expressed  as  a  finite  3-D  tensor  polynomial, 

(50)  a  =  ^C^T^ 

A 


comprised  of  a  linear  combination  of  all  the  independent  tensor  products  (generators)  formed  from  the  three 
primary  tensors.  The  coefficients  in  the  polynomial  are  functions  of  the  independent  invariants  of  the  tensors. 
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For  this  problem  the  dimension  of  the  [minimal]  tensor  base  is  A  =  41  (cf.  Spencer  [43]).  This  is  very  large 
and  is  unlikely  to  be  used  in  practice. 

The  complexity  presented  by  such  a  large  tensor  basis  can  be  side  stepped  by  simplifying  approximation 
regarding  the  63  term.  It  has  been  argued  [21,  27]  that  for  the  range  of  values  used  for  the  constant  C2  the 
inequality  (7$  Q  holds  and  thus  the  term  multiplied  by  will  only  have  a  small  effect  on  the  solution. 
This  approximation  decouples  the  contributions  of  S*  and  R*  to  a.  As  the  equation  is  linear  the  solution 
is  determined  using  the  superposition  principle.  This  one  allows  to  split  the  problem  into  two  equations  of 
lower  tensor  base  dimension: 

(51)  a  = 
where  a^  stands  for  the  solution  dependent  on  S*, 

(52)  a^  =  -QT  [61S*  -  -  fia^)] 

and  a^  denoting  the  solution  dependent  on  R* 

(53)  a^  -  -gr  [62R*  -  64(a^fi  -  fia^)]  . 

The  decomposition  of  a  into  portions  dependent  on  S*  and  R*  is  unique. 

Applying  the  results  from  [21]  the  strain  dependent  portion  of  the  solution  for  the  anisotropy  tensor  can  be 
written 

a®  =  -2C^tS*  -  4a2T^{S*n  -  fiS*)  -  _  ^{S*fi2}5) 

o 

(54)  -16a4r'*(f2S*fi2  -  fj2s*J7)  -  32a5T®{S*n2}(fj2  _ 

o 

where  =  big{l  +  0:2  =  as  =  ^biblg^hi,  04  =  -§61643^/11,  as  =  ^biblg^hi, 

ho  =  b^gr,  hi  =  /i2[l  +  2hoa;^]“^  and  /12  =  [2  +  hQU)'^]~^. 

Using  a  similar  procedure  the  portion  of  the  Reynolds  stress  anisotropy  dependent  on  the  baroclinic  dyad  is 
written 

=  -2C'^rR*  -  4a2r2(R*n  -  OR*)  -  Sasr^C^^R*  +  R*fi^  - 

O 

(55)  -16a4r^(fiR*fi2  -  fl^R*^)  -  32a5T5{R*fi2}(o2  _  i{f^2}5) 

o 

in  which  the  coefficients  have  the  same  form  as  those  in  a*^  with  the  exception  that  bi  is  replaced  by  62.  To 
obtain  the  full  anisotropy  tensor  the  two  complex  expressions  (54)  and  (55)  need  to  be  added.  In  the  light 
of  the  complexity  of  the  three-dimensional  formalism  a  two-dimensional  formalism  is  developed. 

4.2.  A  two-dimensional  algebraic  Reynolds  stress  model.  A  simpler  and  tractable  two-dimensional 
treatment  is  possible.  In  many  engineering  flows  the  mean  flow  and  the  statistics  of  the  turbulence  are  two- 
dimensional.  The  two-dimensional  problem  is  less  complicated  as  the  number  of  tensor  products  necessary 
to  express  the  solution  is  substantially  reduced.  The  symbols  S,  ,  R  now  denote  two-dimensional  tensors. 

It  is  necessary  to  recast  the  equation  for  the  anisotropy  in  terms  of  the  traceless  2-D  tensors:  ^j{{u))  = 

~  J^jipiP)  =  ^j{PiP)  ~  ^^p{PiP)^if  •  Here,  the  two-dimensional  Kronecker  symbol 
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is  <5^^^  =  [^i?]  =  1  for  i  =  j  =  1, 2  and  0  otherwise.  The  pressure-strain  model  is  then  written: 


^ij  -  P  4  =  -C'lP  +  P{^) 
[1  —  C3  +  2^2] 


^ipS.pjii'^))  r^S.pqi{'^))^pq^i3 


hh'" 

2, 

3^ 


SU{u))+ 


[1  —  C4  “  ^d2][(^ipQ.pj{{'^))  ^lip{{'^))^pj]~^ 


-  ^2  Spp  ( (u. ) )  dij 


(56)  2[1- (73  +  2^2] 

The  pressure  acceleration  is  written 

(57) 


Opj  ^pj 


'ip 


(2)’ 

Pj 

2 


■4 

2  1 

Spp{{u)) 

1 

Spg 

3 

J  1 - 1 

1 

■ip\ 

3 

Skk{{u)) 


(2)1 

a 

2 


12  12 

Mij  =  --■^TuMk)Rij  +  ^■^ruUo{k)Rpp 


X  <5^2) 

^ij  ^ij 


The  fact  that  both  2D  and  3D  expressions  of  the  pressure-strain  correlation  model  must  give  the  same  result 
when  applied  to  two-dimensional  mean  flows  will  be  used  for  our  simpliflcations.  Recasting  the  model  in  2D 
is  done  to  take  advantage  of  the  simplifications  that  result  from  the  2D  structure. 

Inserting  the  closures  for  the  pressure-strain  and  mean  acceleration  terms,  (56),  (57)  into  the  Reynolds  stress 
equation,  using  the  same  ansatz  regarding  transport  at  a  fixed  point,  the  2D  analog  of  (48)  is  obtained: 


a  =  -gr 


(58) 


61S*  +  62E*  +  h  ^aS*  +  S*a  -  |{aS*}5)  -  h  (aO  -  fia) 

)-H- 


^^3  2 


d  ,  S 


-  3-2 


with  61  =  62  =  ^|tuI/o,  i>3  =  C3-2d2,  64  =  C4+2d2,  h  =  biSpp{{u))+b2Rpp{p,p),  be  =  2b3Spp{{u)) 

and  g  having  the  same  expression  as  in  the  3-D  algebraic  equation.  The  two-dimensional  polynomial  solution 
of  (58)  is,  as  before,  also  written  as 


(59) 


a  =  53C^T\ 
A 


Unlike  the  three-dimensional  solution,  however,  the  generators  now  consist  of  only  five  tensor  groups: 


(60) 


^ 


5 


T^=S*,  T2  =  S*fi-OS*,  T3  =  R*,  T^  =  R*f2-fiR* 


for  which  there  are  five  non-zero  independent  invariants, 

(61)  a2  =  {S*^},  a;2  =  -{0^},  {R*^},  {R*S*},  {S*R*S1}. 

The  fact  that  there  are  no  other  independent  tensor  generators  [or  invariants]  can  be  verified  using  the  2x2 
matrix  identity: 


(62) 


2abc  =  bc{a}  +  a{bc}  +  ac{b}  -  b{ac}  +  ab{c}  +  c{ab}  — 
c{a}{b}  -  a{b}{c}  +  ({ac}{b}  -  {acb})5^^\ 
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To  obtain  the  solution  to  the  algebraic  equations  for  the  anisotropy  tensor,  (58),  a  procedure  similar  to  the 
one  devised  by  [20]  is  used.  Three  5x5  matrices  are  defined: 


(63) 


ijirj 


S  _ 

3  2  J 


T^Q*  -  n*T^  =  ^ 

A 


A 


A 


for  which  A  =  0  —  4.  The  elements  of  the  matrices  are  determined  from  the  above  equations  by  making  use 
of  matrix  relations  stemming  from  the  Cayley-Hamilton  theorem.  The  2-D  tensor  polynomial  a  =:  XIa 
is  introduced  in  both  sides  of  relation  (58).  By  making  use  of  the  above  matrix  identities  the  coefBcients  of 
the  tensor  polynomials  are  found  to  satisfy  the  following  system  of  equations: 


(64) 


=  -gr 


+  ^»2^3A  +  f>3  ^ cmi - 64 5] - bsSox -hY, 


The  solution  of  this  system  of  equations  determines  the  model  coefficients  in  the  following  algebraic  expression 
for  the  anisotropy  tensor 

5  -  d(2))  +  Q2h9fir{s*n  -  n  s*) 

(65)  -2ClT[R*  +  h9fiT{R*n-aR*)] 

suitable  for  two-dimensional  mean  flows.  The  eddy  viscosities,  and  are  given  by: 


a  =  —2C^t 


Q2S*  -f  (Qi  -f  Q3)b3gf2Ta^ 


(66) 

(67)  c:,  = 

The  Qi  coefficients  are  given  by 

(68) 

(69) 


big  f  1/2 


1  -  |(i>3ffT)2(72/i/2  +  2(645/iT)2a;2 


b2gfi/2 


l-f  2(643/1  r)2w2  ■ 

{S*R*} 


61 


*T,*i  {s*R*m 

— -h2g/iT64^ 


O  1  I  ^  (^3gT)VVi/2  ,  .N  ^3^5  . 

3 1-1-2(643/1  t)2w2^^^  ^  3bl^^^ 


(70) 

(71) 


65  1 -f  2(643/iT)2a;2 

bih  2p/ir<T2 


Q3  =  - 


with  /i  =  (1  +  begr/d)  ^  and  /2  =  (1  —  bQgT/6)  Note  that  the  direct  effect  of  compressibility  as  reflected 
in  the  baroclinic  dyad  occurs  in  Qi  and  Q2. 


The  high  order  of  nonlinearity  of  the  algebraic  equations  does  not  permit,  in  general,  the  construction  of  a 
fully  explicit  solution.  Instead,  an  iterative  approach  is  employed  during  the  computations  to  generate  the 
correct  values.  The  algebraic  solution  is  linearized  by  lagging  the  turbulence  production  term  which  contains 
the  nonlinearity. 
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4.3.  Discussion.  To  conclude  this  section  some  general  statements  regarding  the  behavior  of  the  al¬ 
gebraic  closure  derived  for  the  Reynolds  stresses  are  highlighted.  In  §3  a  closure  for  the  Reynolds  stress 
equations  was  obtained.  In  the  present  section  the  fixed  point  solution  of  the  modeled  second  moment 
equations,  under  the  condition  of  structural  equilibrium,  was  obtained.  For  two-dimensional  mean  flows  the 
compressible  algebraic  stress  model  can  be  symbolically  written  as 

(72)  a  =  -  I/to  r  -  *^11  [r  n  i/t2[^5  -  5^"^]  -  I'^o  S*  “  -  «  E*]- 


Note  that  products  of  the  mean  strain  and  the  baroclinic  dyad  do  not  appear.  This  is  due,  for  two-dimensional 
flows,  to  the  relation  T^S*  -b  -  |{T'?S*}(5  =  -2{T’?S*}T°;  the  generator  comprised  of  the  mean  strain 
and  baroclinic  dyad  product  is  redundant.  Examining  the  above  expression,  the  following  observations  can 
be  made: 


1.  The  first  two  terms  are  the  same  terms  obtained  in  algebraic  stress  closure  for  two-dimensional  [in 
the  mean]  incompressible  flows. 

2.  The  first  two  eddy  coefficients  are  functions  of  the  relative  strain  and  the  relative  rotation  as  is  the 
case  in  incompressible  flows.  They  are  now  also  functions  of  the  turbulent  Mach  number  and  the 
gradient  Mach  number. 

3.  Neither  of  these  eddy  coefficients  depend  on  baroclinic  effects. 

4.  It  is  seen  that  in  the  absence  of  mean  velocity  gradients  that  the  turbulence  is  anisotropic  due  to  the 
mean  baroclinic  dyad.  This  anisotropy  manifests  itself  in  the  deviatoric  as  well  as  diagonal  terms. 


It  is  useful  to  construct  a  simple  example  to  see  how  the  new  effects  influence  the  anisotropy.  Consider  a 
simple  shear  flow  with  a  streamwise  acceleration:  let  the  Favre  mean  velocity,  the  mean  density  and  the 
mean  pressure  gradients  be  represented  by  {71,2,  i7i,i,  Vp  =  /9,2,  and  VP  =  P,i. 


(73) 

(74) 

(75) 

(76) 


2  11 

2  ^  ^152  P»2  P»1 

111 

022=  2^*1  ^1,2  +  ^^t2  —  ^tiUi^2  P,2Pa 

O33  =  2  “  3  ^*2 

012  =  -r  2>'to!7i,2  !7i,i  Pjl 


Several  observations  regarding  the  above  algebraic  expression  for  the  anisotropy  can  be  made: 


1.  The  expression  is  the  first  [that  we  know  of]  rigorous  indication  of  the  direct  role  the  baroclinic  dyad 
plays  in  determining  the  Reynolds  stresses. 

2.  The  expression  indicates  that,  for  arbitrary  mean  deformation,  the  mean  baroclinic  dyad  contributes 
to  the  deviatoric  portions  of  the  Reynolds  stress. 

3.  The  expression  also  indicates  that  the  baroclinic  dyad  also  changes  the  relative  magnitude  of  the 
normal  stresses.  This  effect  only  occurs  for  mean  deformations  that  are  rotational.  For  an  irrotational 
mean  deformation  the  baroclinic  dyad  makes  no  contributions  to  the  normal  stresses. 

4.  For  a  uniform  mean  velocity  the  baroclinic  dyad  is  a  source  of  anisotropy  but  only  in  the  deviatoric 
portions  of  the  anisotropy. 

5.  The  expression  indicates  the  inapplicability  of  any  heuristic  gradient  transfer  arguments  for  the 
Reynolds  stresses  in  flows  with  important  gradients  of  mean  density  and  pressure. 
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While  these  results  indicate  the  inapplicability  of  any  form  of  eddy  viscosity  model  for  the  stresses  in 
compressible  flows  with  arbitrary  large  density  and  pressure  gradients  some  qualifications  are  in  order.  The 
presence  of  the  baroclinic  dyad  is  likely  to  be  important  only  in  rapidly  accelerating  aerodynamic  flows  or 
in  combusting  flows  where  one  can  expect  the  mass  fluxes  to  be  important.  In  the  absence  of  these  effects  it 
appears  that  the  parameterization  of  the  Reynolds  stresses  in  terms  of  powers  of  the  mean  deformation  with 
modifications  according  to  the  compressibility  of  the  fiuctuations  as  those  indicated  in  §3  is  appropriate. 

5.  Computational  investigations  of  free  shear  flows.  The  theory  and  results  presented  in  the 
previous  sections  are  now  implemented  over  a  very  wide  range  of  mean  flow  Mach  numbers  in  the  simulation 
of  free  shear  flows.  Simulations  using  second-order  moment  (SOM)  closures  as  well  as  the  quasi-explicit 
algebraic  models  are  conducted.  The  numerical  experiments  are  constructed  with  the  intention  of  investi¬ 
gating  several  different  issues  of  relevance  to  the  prediction  of  compressible  turbulent  flows  for  engineering 
purposes. 

In  addition  to  assessing  the  sensitivity  of  the  pressure-strain  model  to  the  different  pressure-dilatation  models 
we  compare  the  computational  results  to  what  is  expected  from  laboratory  and  numerical  experiments.  Of 
particular  interest  is  the  well  known  effect  of  compressibility  on  reduction  of  the  spread  rate  of  the  mixing 
layer.  The  main  objective  is  to  assess  the  effectiveness  of  the  compressibility  corrections  in  reproducing 
the  reduced  growth  rate  of  the  free-shear  layers  with  increasing  Mach  number  -  a  phenomenon  which  is 
well  documented  experimentally  [44].  We  also  study  the  adequacy  of  the  algebraic  stress  closure  with  the 
results  of  the  SOM  simulations.  In  particular  the  Reynolds  stresses  and  the  anisotropy  computed  using  both 
of  these  procedures  are  compared  to  solutions  without  the  compressible  corrections  in  order  to  assess  the 
effects  of  the  new  pressure-strain  closure  on  the  anisotropy.  Our  intention  is  to  further  investigate  the  notion 
that  the  dramatic  reduction  of  the  mixing  layer  growth  rate  is  due  to  the  effects  of  compressibility  on  the 
pressure-strain  and,  consequently,  its  effect  on  the  reduction  of  the  turbulence  shear  stress. 

5.1.  Numerical  method.  The  simulations  are  conducted  using  a  finite  difference  scheme,  following 
[45] ,  second-order  accurate  in  time  and  fourth-order  accurate  in  space.  The  governing  equations  are  integrated 
explicitly  in  time  using  the  predictor-corrector  finite  difference  scheme  developed  by  Gottlieb  and  Turkel  [46]. 
The  Gottlieb-Turkel  scheme  is  a  higher  order  accurate  variant  of  the  MacCormack  [47]  predictor-corrector 
method.  During  a  numerical  sweep,  the  inviscid  fluxes  are  alternatively  differenced  backward  in  the  predictor 
step  and  forward  in  the  corrector  step.  Fourth-order  central  differences  are  used  for  the  viscous  and  heat  flux 
terms  as  well  as  for  the  derivatives  in  the  source  vector.  To  maintain  stability,  it  is  required  that  the  CFL 
number  be  less  than  2/3.  To  prevent  numerical  oscillations  in  regions  of  large  gradients  a  smoothing  scheme 
devised  by  MacCormack  and  Baldwin  [48]  is  employed.  The  method,  as  outlined  in  [45],  adds  artificial 
viscosity  that  is  very  small  everywhere  except  in  the  regions  where  the  pressure  oscillates. 

All  simulations  are  conducted  on  a  uniformly  spaced  grid  in  the  computational  domain.  By  means  of 
a  coordinate  transformation  the  mesh  is  transversally  compressed  in  the  physical  domain  in  the  region 
corresponding  to  the  mixing  layer. 

The  physical  domain  is  a  rectangular  box  defined  by  the  set  of  points  {x^y)^  in  which  x  represents  the 
streamwise  coordinate  and  y  the  transversal  coordinate.  The  grid  overlaying  the  computational  domain  of 
size  x5u)  X  ydijj  has  100  x  60  points,  where  the  vorticity  thickness  —  {u\  —U2)/{d{u)/dy)max^  the  brackets 
0  denoting  the  averaged  value.  In  this  nondimensionalization,  the  reference  length  scale  is  the  magnitude 
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of  the  vorticity  thickness  at  the  inlet.  The  initial  profile  for  the  mean  axial  velocity  is  adjusted  such  that 
the  resulting  nondimensional  vorticity  thickness  at  the  inlet  is  equal  to  one.  The  values  of  the  physical 
dimensions  are  {x,  y)  =  [80, 20]  for  the  mixing  layer.  To  evaluate  the  grid  sensitivity  the  number  of  grid 
points  is  increased  by  a  factor  of  1.5.  The  change  in  the  steady-state  values  of  the  peak  of  the  turbulent 
stresses  is  less  than  2  percent. 

To  accelerate  convergence  to  the  stationary  state,  a  local  time  stepping  technique  is  used.  The  convergence 
criterion  is  imposed  so  that  the  global  averaged  residual  profile  attained  is  stationary  for  each  dependent 
variable.  Although  more  stringent  criteria  can  be  sought,  it  is  known,  [49],  that  predictor-corrector  schemes 
are  limited  in  their  ability  to  achieve  very  high  rates  of  residual  reduction. 

Due  to  the  nonlinearity  of  the  destruction  terms,  the  k  -  e  equations  display  a  stiffness  which  can  either 
generate  numerical  instabilities  or  increase  the  computational  time.  In  order  to  avoid  these  inconveniences 
the  turbulence  source  terms  are  treated  implicitly.  The  k  —  e  destruction  terms  are  decoupled  by  suitable 
manipulation  of  the  e/k  ratio,  considered  as  a  known  quantity  from  the  previous  time  step. 

Initial  and  boundary  conditions.  The  initial  fields  are  obtained  by  propagating  the  inflow  conditions 
throughout  the  entire  domain;  hence,  the  flow  has  to  sweep  the  domain  at  least  one  time  to  obtain  meaningful 
results.  For  ease  of  computation,  the  inflow  initial  conditions  (IC)  for  the  flow  variables  are  assumed  to  be  a 
smoothed  step  (hyperbolic  tangents  for  the  axial  velocity  or  species)  or  a  smoothed  hat  profile  (for  turbulence 
quantities).  Uniform  profiles  are  assigned  at  the  inlet  to  the  static  pressure,  to  the  static  temperature  and, 
by  virtue  of  the  equation  of  state,  to  the  static  density. 

The  boundary  conditions  (BC)  are  set  according  to  the  elliptic  nature  of  the  problem  on  all  four  boundaries. 
The  inflow  BC  are  fixed  (Dirichlet)  for  all  primary  variables  in  the  supersonic  and  subsonic  regions  with  one 
exception.  For  the  portions  where  the  flow  is  subsonic  the  pressure  is  allowed  to  adjust  to  the  characteristic 
waves  through  a  Neumann  boundary  condition.  At  the  outflow  and  outer  boundaries  zero  gradient  (Neu¬ 
mann)  conditions  are  applied  due  to  their  non-reflective  properties  in  relation  with  the  outgoing  waves.  In 
the  mixing  layer,  the  static  pressure,  temperature  and  density  in  the  two  free  stream  layers  are  matched. 
This  isolates  the  contributions  to  the  layer  growth  to  those  due  solely  to  the  variation  of  the  reference  Mach 
number. 

Pree-shear  flow  parameters.  The  magnitude  of  the  effects  of  compressibility  are  often  parameterized  by 
the  convective  Mach  number  [50].  In  our  formulation  the  expression  for  the  convective  Mach  number  is 
Me  =  M(1  -  ry)/2,  where  =  U2/U1  is  the  axial  velocity  {u)  ratio.  A  wide  range  of  values  of  Me  including 
both  subsonic  and  supersonic  regimes  are  considered:  0.2  <  Me  <  2.0.  The  convective  Mach  number  is 
varied  by  keeping  the  velocity  ratio  constant  and  varying  the  reference  Mach  number  M.  The  reference 
Reynolds  number  is  Re  =  b  x  10®,  while  all  the  other  nondimensional  parameters  are  kept  constant  well. 
The  subscript  1  refers  to  the  high-speed  stream  value,  while  the  subscript  2  refers  to  the  low-speed  stream 
value  in  the  mixing  layer.  No  jet  simulations  are  shown  in  this  article;  they  can  be  found  in  [37]. 

In  the  results  presented,  the  spatial  coordinate,  for  the  hydrodynamic  variables,  are  given  by  r).  In  the 
mixing  layer,  rj  —  ^  where  u  =  -^73™-  Unless  otherwise  specified,  the  results  shown  are  from 

computations  [ARSM  and  SOM]  with  models  based  on  the  dilatational  models  of  [8]. 

After  an  initial  period  of  flow  development  a  linear  growth  rate  of  the  shear  layers  is  attained.  In  the  fully 
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developed  regime,  when  a  linear  growth  rate  has  been  established,  the  mean  velocity  and  the  normalized 
Reynolds  stresses  display  self-similar  behavior.  The  spread  rate  of  a  turbulent  shear  layer,  in  the  self-similar 
region,  is  conventionally  expressed  as  dSu/dx  =  Cs{l  -  r^) / {1  +  r^)  where  5u  could  be  either  the  10  percent 
visual  thickness  of  the  shear  layer  based  on  the  normalized  velocity  profile  or  the  vorticity  thickness,  and 
Cs  is  a  constant  (approximately).  As  the  present  article  is  concerned  with  compressible  corrections  to  the 
turbulence  moment  equations  no  plots  of  the  mean  flow  quantities  are  given;  the  mean  field  variables  are 
quite  typical  and  can  be  seen  in  the  literature,  [11,  51,  52]. 

5.2.  Mixing  layer  simulations.  The  primary  concerns  of  this  article  are  with  1)  a  representation  of 
the  effects  of  compressibility  as  they  appear  in  the  pressure-strain  covariance  and  2)  the  construction  of  an 
algebraic  Reynolds  stress  model  useful  for  engineering  calculations.  Studies  related  to  both  these  issues,  for 
the  Me  =  1.07  mixing  layer  configuration,  are  shown  in  the  first  two  figures.  The  Reynolds  stresses  and 
the  production/dissipation  ratio  obtained  for  ARSM  and  SOM  calculations  are  shown  in  Figures  1  and  2. 
Also  shown,  as  a  double  check,  are  curves  based  on  an  a  priori  evaluation  of  the  algebraic  stress  model; 
the  data  from  a  SOM  calculation  is  used  as  input  to  the  ARSM.  The  neglect  of  turbulent  transport  in  the 
algebraic  stress  model  is  reflected  in  differences  in  the  normal  stresses  at  the  centerline.  In  general  the  a 
priori  evaluated  ARSM  curves  follow  closely  the  computed  ARSM  values  indicating  that  for  this  flow  the 
neglect  of  transport  in  the  moment  transport  equations  is  justified.  The  algebraic  approximation  is,  as  has 
been  verified  at  several  different  convective  Mach  numbers,  suitable  for  an  investigation  of  trends  in  this  flow 
configuration. 

Also  shown  in  Figures  1  and  2  is  the  SOM  computation  without  compressibility  corrections  [NC  -  no  com¬ 
pressibility].  The  effect  of  the  compressibility  corrections  manifest  themselves  as  decreases  in  the  centerline 
values  of  {uu)^  (uv)  and  (vv)  of  approximately  16%,  25%  and  24%  respectively.  This  is  consistent  with  that 
seen  in  DNS.  In  addition,  the  streamwise  normal  stress  suffers  a  smaller  reduction  than  the  other  Reynolds 
stresses.  This,  as  will  be  seen,  will  manifest  itself  in  an  increase  in  the  streamwise  normal  anisotropy. 

The  production/dissipation  ratio  is  shown  in  the  bottom  of  Figure  2.  The  average  level  of  P/e  for  the 
compressible  case  is  smaller  than  the  incompressible  case;  this  is  consistent  with  the  reduced  turbulent 
kinetic  energy  of  the  compressible  flow.  While  the  overall  level  of  the  production  is  smaller  it  is  interesting 
to  note  that  at  one  location  in  the  mixing  layer,  the  centerline.  P/e  is  nominally  higher.  It  has  been 
conjectured  that  this  might  be  due  to  the  higher  relative  strain  rates  compressible  flows  may  be  able  to 
sustain.  This  is  an  issue  that  further  DNS  may  resolve. 

Figure  3  displays  the  vorticity  thickness.  Vorticity  thickness  spread  rates  as  predicted  by  the  SOM  simula¬ 
tions,  with  and  without  (NC)  the  compressibility  correction,  are  shown  on  the  top  of  Figure  3.  It  is  seen, 
that  even  with  such  modest  contributions  from  compressibility,  scaling  as  they  do  with  a  sizable  change 
in  growth  rate  is  seen.  The  variation  of  Mt  with  Me  is  shown  in  Figure  5. 

Figure  3  also  displays  the  different  components  of  the  pressure-strain  tensor  computed  from  the  SOM  simula¬ 
tion  with  and  without  (NC)  the  compressibility  correction.  The  different  components  of  the  pressure-strain 
are  all  reduced  by  about  the  same  amount.  This  is  precisely  the  behavior  seen,  in  both  direction  and 
magnitude,  in  the  recent  DNS  [15]  -  see  Figure  4,18. 

Reynolds  stress  anisotropy.  In  Figure  4  the  computed  values  of  the  centerline  Reynolds  stress  anisotropies 
are  shown.  At  low  values  of  compressibility  the  anisotropy  has  values  similar  to  that  of  simple  incompressible 
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shears.  At  higher  levels  the  anisotropy  behaves  very  much  like  the  compressible  shear  DNS;  more  of  the 
energy  of  the  turbulence  is  in  the  streamwise  component  as  made  clear  by  the  increase  in  an.  This  is  at  the 
expense  of  the  crosstream  component,  (vv),  which  is  reduced.  This  is,  as  shall  be  discussed  further,  due  to 
the  reduction  of  intercomponent  energy  transfer  due  to  the  compressible  aspect  of  the  pressure-strain.  There 
is  also  a  reduction  in  the  shear  stress  as  is  seen  in  compressible  DNS. 

Comparisons  with  laboratory  Reynolds  stresses.  The  effects  of  the  new  compressible  models  on  the 
anisotropy  is  shown  in  Figure  4.  Unfortunately,  in  the  laboratory  situation,  the  spanwise  Reynolds  stresses 
are  rarely  measured.  As  a  consequence  the  anisotropy  is  not  known  and  one  returns  to  primitive  variables. 

The  maximum  values  of  diverse  Reynolds  stresses  in  the  mixing  layer  are  plotted  as  function  of  the  convective 
Mach  number  in  Figure  4.  In  this  figure  and  ay  represent,  respectively,  the  (uu)  and  (vv)  Reynolds  stresses 
nondimensionalized  by  the  square  of  the  velocity  difference  across  the  layer.  Also  shown  are  the  experimental 
results  of  Elliot  and  Samimy  [53].  Underlying  this  comparison  is,  of  course,  the  implicit  assumption  of 
self-similarity.  The  maximum  Reynolds  stresses  decrease  with  compressibility.  In  the  present  simulations, 
however,  the  peak  values  of  the  turbulence  intensities  (axial  and  transverse)  do  not  vary  as  strongly  as  in 
the  laboratory  findings  of  Goebel  and  Dutton  [54].  In  fact  the  present  simulations  are  closer  to  the  results 
of  either  [54]  or  [53]  than  either  [54]  or  [53]  are  to  each  other. 

There  appears  to  be  a  discrepancy.  A  survey  of  the  diverse  compressible  DNS  indicates  that  an  increases  with 
compressibility.  This  implies  that  {uu)  decreases  faster  than  (vv).  [Figure  9  of  [14]  is  the  most  comprehensive 
presentation  of  this  trend.]  This  does  not  appear  to  be  the  case  in  the  laboratory  experiments.  In  the  mixing 
layer  of  [53]  (uu)  decreases  more  rapidly  than  (vv)  for  all  Me  implying  that  an  decreases  with  increasing 
compressibility  contrary  to  DNS  results.  In  [54]  (uu)  decreases  more  rapidly  than  (vv)  for  low  Me]  while 
at  moderate  Me,  (uu)  decrease  less  rapidly  than  {vv)  in  accord  with  DNS  results.  No  explanation  for  this 
discrepancy  between  laboratory  and  numerical  data  is  known. 

The  pressure-strain  tensor.  The  effects  of  compressibility  on  the  Reynolds  stresses  are  believed  to  be  the 
source  of  the  dramatic  influence  on  the  spread  rate  as  the  convective  Mach  number  increases.  The  physical 
mechanisms  responsible  for  these  phenomena  were  once  believed  to  be  the  compressible  dissipation  and 
pressure-dilatation  which  caused  a  decrease  in  the  Reynolds  stresses  due  to  an  overall  decrease  in  k.  This 
was  accompanied  by  no  substantial  changes  in  the  anisotropy.  As  discussed  in  §1,  recent  simulations  suggest 
that  compressibility  causes  a  turbulent  shear  stress  reduction  due  to  a  reduction  in  its  associated  anisotropy. 
This  is  a  statement  that  compressibility  manifests  itself  structurally,  not  energetically.  The  reduction  in  the 
anisotropy  leads,  of  course,  to  a  reduction  in  k  and  ultimately  in  the  mixing  layer  growth  rate. 

More  recent  simulations,  [9,  14,  15,  16]  all  seem  to  indicate  that  the  major  changes  in  the  shear  stress  are 
due  to  the  effect  of  compressibility  on  the  pressure-strain  covariance.  Figure  6  is  an  illustration  of  this  idea: 
Figure  6  depicts  the  different  components  of  pressure-strain  tensor  and  its  trace  [the  pressure  dilatation]  at 
the  centerline  of  the  mixing  layer.  Figure  6  is  patterned  after  Vreman  et  al  ^s  [14]  Figures  3  and  9.  The 
behavior  of  all  the  components  of  the  pressure-strain,  shown  in  Figure  6,  are  very  much  in  accord  with  the 
figures  of  [14].  The  only  notable  difference  is  the  relative  levels  of  1122  and  lisa.  The  trends  are  very  much 
the  same  and  the  relative  behavior  -  a  commensurate  across  the  board  reduction  with  Me  -  is  very  much  in 
accord  with  [15]. 
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The  behavior  of  the  pressure-strain  illustrated  in  Figure  6  is  very  much  in  keeping  with  the  current  under¬ 
standing  of  the  mechanism  for  the  reduction  of  the  turbulence  energy  by  the  reduction  of  the  shear  stress: 
the  only  source  for  the  kinetic  energy  is  in  the  {uu)  component  and  its  production  is  proportional  to  {uv). 
The  reduction  [in  magnitude]  of  Hu  and  1122  means  that  less  energy  generated  in  the  (uu)  equation  is 
transferred  to  {vv)  [by  1122]  ■  As  (vv)  has  no  production  (vv)  is  much  smaller.  As  consequence  of  the  fact 
that  the  production  of  {uv)  is  proportional  (vv)  much  less  (uv)  is  produced  and  much  less  (uu)  is  produced. 
This  can  be  seen  by  the  considering  the  truncated  forms  of  the  second-order  moment  equations: 


(77) 

^  {uu)  ~  -  {uv)Ul,2  +  IIii  + 

(78) 

^  {vv)  ~  1122  +  •••• 

(79) 

^  {uv)  ~  -  {vv)Ul,2  -1-  ni2  -1- 

The  thesis  of  [15]  provides  a  schematic  of  the  above  process.  The  behavior  of  the  Reynolds  stresses  and 
the  pressure-strain  indicated  by  Figures  2  and  6  is  consistent  with  the  mechanism  for  the  reduction  just 
delineated.  Specifically,  the  fact  that  less  of  the  {uu)  energy  is  transferred  to  the  other  components  of  the 
Reynolds  stress  means  that  the  anisotropy,  an  =  26ii,  must  increase  as  seen  in  all  the  DNS.  Note  also 
that  the  pressure-dilatation  as  derived  by  the  asymptotic  procedure  in  [8]  is  extremely  small;  negligible 
in  comparison  to  the  pressure-strain.  Yet  its  use  to  obtain  the  deviatoric  portions  of  the  pressure-strain 
produces  a  very  sizable  reduction  in  the  shear  stress  and  the  growth  rate. 

Gradient  Mach  number.  Sarkar  [16]  and  earlier  work  cited  therein  have  drawn  attention  to  the  gradient 
Mach  number  as  a  potentially  useful  parameterization  of  the  effects  of  compressibility.  Many  of  the  effects  of 
compressibility,  as  indicated  by  diverse  DNS,  have  been  observed  to  become  more  apparent  as  the  gradient 
Mach  number  increases.  Shown  in  the  bottom  of  Figure  5  is  the  gradient  Mach  number  using  the  definition 
Mg  =  |Mta(fc)/es.  The  definition  of  the  gradient  Mach  reflects  the  fact  that  the  Kolmogorov  scaling, 
e  ~  {k)^^^/£  has  been  used  to  eliminate  the  length  scale  in  the  definition  Mg  =  aijc.  In  DNS,  the  length 
scale  i  is  determined  from  the  two-point  correlation;  no  such  opportunity  occurs  in  single-point  closures. 
Both  the  DNS,  for  example  Figure  4.29  in  [15],  and  the  SOM  simulation  indicate  a  similar  decrease  in  rate 
of  increase  of  Mg  with  Me-  The  decrease  seems  to  occur  at  lower  Me  in  the  DNS.  This  is  likely  to  be  related 
to  the  different  length  scales  used. 

One  of  the  noteworthy  observations  made  of  the  effects  of  compressibility  is  their  tendency  to  saturate,  [11]. 
In  the  bottom  of  Figure  5  the  rate  of  increase  of  the  maximum  turbulent  Mach  decreases  at  higher  convective 
Mach  numbers.  This  is  line  with  observations  of  [11]. 

Kolmogorov  scaling.  While  in  DNS  a  length  can  be  determined  firom  the  two-point  correlation,  no  such 
possibility  exists  for  single-point  closures.  Yet  a  length  scale,  reflecting  the  two-point  nature  of  the  turbulence 
problem,  is  required.  The  length  scale  appearing  in  the  Kolmogorov  inertial  range  scaUng,  e  ~  {k)^^‘^/£,  was 
used  in  [8]. 

There  is  a  proportionality  coefficient,  a,  in  the  Kolmogorov  scaling  a  =  est/{2k/2)^/^.  For  infinite  Reynolds 
number  isotropic  turbulence  q  ~  1.  For  finite  but  large  Reynolds  number  anisotropic  turbulence  undergoing 
deformation  the  Kolmogorov  scaling  is  likely  to  be  useful  but  a  is  likely  to  be  a  flow  dependent  quantity. 
Sreenivasan  [55]  has  made  some  studies  of  the  variation  of  the  Kolmogorov  scaling  coefficient:  [55]  has  found 
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values  of  a  ranging  between  0.4  -  2  for  diverse  incompressible  simple  shear  flows.  To  allow  for  the  expected 
variability  of  a  two  different  values  of  are  used  in  many  of  the  simulations.  The  value  of  a  does  not  change 
any  of  the  trends  and  does  not  seem  to  have  too  strong  an  effect  on  the  energy  at  the  centerline  as  seen  in 
Figures  4  and  5.  It  does  however  effect  the  mixing  layer  spread  rate  as  can  be  seen  in  Figure  8. 

Mixing  layer  growth  rates.  Settles  and  Dodson  [44]  have  compiled  a  very  large  number  of  experimental 
results.  These  are  shown  in  Figures  7  and  8.  The  experimental  curves  come  from  the  experiments  of 
[50,  54,  56,  57,  58,  59,  60,  61].  The  scatter  in  the  data  reflects  the  fact  that  different  experiments  are  done  in 
different  wind  tunnels,  with  different  inlet  configurations  and  with  different  reservoirs.  The  computational 
curves  shown  in  Figures  7  and  8  reflect  two  different  but  complementary  investigations.  As  the  “Langley 
curve”  [62]  has  been  a  popular  benchmark  it  has  also  been  included.  The  well  known  reduction  in  the  mixing 
rate  with  increasing  Me  is  seen  in  all  the  data. 

Figure  7  reflects  computations  using  the  “energetic”  approaches.  A  k  —  e  scheme  has  been  used  to  calculate 
the  mixing  layer  growth  and  the  only  compressible  corrections  are  in  the  k  equation;  the  Reynolds  stresses  are 
closed  with  the  usual  incompressible  eddy  viscosity  [Boussinesq]  approximation.  The  compressible  dissipation 
and  pressure-dilatation  models  of  Zeman  [4],  Sarkar  et  al  [6,  38],  Ristorcelli  [8]  are  used  in  the  fc-e  scheme. 
Figure  7  indicates  that  the  results  based  on  the  models  of  Zeman  [4],  and  Sarkar  et  al  [6,  38]  capture 
the  mixing  layer  reduction.  The  reduction  of  the  mixing  layer  growth  rate  has  been  captured  by  similar 
mechanisms  in  both  models;  most  of  the  suppression  is  provided  by  a  substantial  amount  of  additional 
dissipation  that  comes  from  the  compressible  dissipation  models.  As  was  pointed  out  in  §1,  the  compressible 
dissipation,  as  has  been  indicated  by  DNS  and  by  asymptotic  analysis,  is  not  important  in  these  flows.  Any 
arbitrarily  dissipative  term  added  to  the  k  equation,  appropriately  calibrated  and  scaled,  would  produce  such 
agreement  for  the  mixing  layer  growth.  Moreover,  as  with  all  energetic  approaches,  there  is  no  possibility  of 
accounting  for  the  very  important  structural  changes  that  appear  in  the  anisotropy. 

As  indicated  by  Figure  7  the  pressure-dilatation  model  of  [8],  with  a  =  2,  accounts  for  a  nominal  suppression 
of  the  growth  rate.  The  compressible  dissipation  being  negligible  in  analysis  of  [8]  is  not  included  in  this 
calculation.  From  incompressible  DNS,  [63]  it  can  be  argued  that  a  value  of  a  «  1  might  be  more  appropriate. 
As  the  pressure  dilatation  scales  with  a  value  oi  a  would  decrease  the  pressure-dilatation  effects  by 
a  factor  of  four. 

Figure  8  reflects  a  computation  using  the  compressible  algebraic  Reynolds  stress  model.  This  is  exactly  the 
same  computation  as  given  in  the  Figure  7  with  the  exception  that  the  algebraic  Reynolds  stress  approxi¬ 
mation  now  includes  the  effects  of  compressibility  in  the  pressure-strain  covariance  representation  used  to 
close  the  second-moment  equations.  Note  the  substantial  changes  on  the  layer  growth  rate  prediction  that 
include  the  effects  of  compressibility  on  the  Reynolds  stresses  produces.  The  change  is  more  drastic  for  the 
Ristorcelli  modeling  than  the  Sarkar  modeling;  this  is  to  be  expected  as  the  majority  of  the  growth  rate 
reduction  using  the  earlier  Sarkar  modeling  was  built  into  a  dissipative  term  and  the  relative  change  is  small. 

The  algebraic  model  built  upon  dilatational  closures  of  [8]  shows  improvement  over  the  results  from  a  simple 
k  —  e  scheme.  Figure  8  depicts  the  outcomes  of  the  computations  with  the  Ristorcelli  [8]  based  algebraic 
closures  for  two  values  of  the  Kolmogorov  scaling  coefficient.  The  results  for  a  =  0.4  with  both  the  pressure- 
strain  and  pressure-dilatation  provide  a  modest  decrease  in  the  spread  rate;  handily  matching  that  predicted 
with  simply  the  pressure-dilatation  with  a —  2.  The  ARSM  calculation  is  substantially  changed  when  a  =  2; 
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if  Oi  were  an  undefined  constant  one  might  be  tempted  to  set  it  so  as  to  match  the  reduction  in  spread  rate 
indicated  by  the  data.  Given  that  little  is  known  about  a  [though  it  is  sure  to  be  in  the  range  0.4-2.0  and 
most  likely  a  ss  1]  these  issues  must  be  left  to  future  DNS  work. 

It  is  concluded  that  the  present  pressure-strain  modeling  method  based  on  the  extension  of  the  well- 
established  incompressible  procedure,  in  which  the  trace-free  constraint  is  relaxed,  does  rationally  account 
for  a  significant  portion  of  the  reduction  in  shear  stress  and  growth  rate.  This  procedure  is  considered  a 
leading  order  contribution  to  the  compressible  turbulence  shear  stress  problem.  The  possibility  that  addi¬ 
tional  compressible  corrections  to  the  pressure-strain,  addressing  physics  not  able  to  be  accounted  for  using 
this  constitutive  relation  methodology,  may  require  consideration. 

6.  Some  issues  in  compressible  turbulence  modeling.  In  our  effort  to  obtain  a  closure  for  the 
effects  of  compressibility  a  few  issues,  already  alluded  to,  have  become  clearer.  These  issues  are  now  high¬ 
lighted. 

The  baroclinic  dyad.  The  procedure  invoked  to  obtain  the  algebraic  Reynolds  model  indicates  that  the 
baroclinic  dyad  is  a  source  of  turbulence.  This,  of  course,  can  be  anticipated  by  inspection  of  the  second- 
order  moment  equations.  The  consequence  is  that  an  eddy  viscosity  representation  for  the  Reynolds  stresses 
is,  from  first  principles,  inappropriate  for  classes  of  flows  in  which  the  mass  fluxes  are  important. 

The  presence  of  the  baroclinic  dyad  is  likely  to  be  important  only  in  rapidly  accelerating  aerodynamic  flows 
such  as  through  shocks  or  in  hypersonic  situations.  The  baroclinic  dyad  is  also  likely  to  be  important  in 
flows  where  one  can  expect  the  mass  fluxes  to  be  important.  In  noncombusting  supersonic  flows  it 
appears  that  a  parameterization  of  the  Reynolds  stresses  in  terms  of  the  mean  deformation  with  modifications 
according  to  the  compressibility  of  the  fluctuations  as  those  derived  in  §3.1  appears  appropriate. 

A  length  scale  for  single-point  models  of  compressible  turbulence.  Much  of  the  work  described  here 
is  based  on  the  leading  order  asymptotic  analysis  of  [8].  In  [8]  the  effects  of  the  non-zero  divergence  of 
compressible  turbulence  was  parameterized  by  several  two-point  integrals  made  nondimensional  by  a  length 
scale.  For  single-point  turbulence  closures  the  length  scale  is  typically  determined  by  the  Kolmogorov  inertial 
range  scaling,  e  ~ 

For  compressible  turbulence  of  interest  to  supersonic  aerodynamic  flows  the  cascade  mechanism  will  be 
comprised  of  the  usual  nonlinear  solenoidal  modal  interactions.  The  Kolmogorov  scaling  is  expected  to 
be  valid  in  the  weakly  compressible  limit.  However,  given  the  notable  behavior  seen  in  simulations  with 
increasing  compressibility  of  a  length  scale,  see  Figure  4.28  of  the  recent  [15],  one  might  expect  the  coefficient 
of  proportionality,  what  we  have  called  a,  to  be  a  function  of  compressibility.  That  this  might  be  the  case 
is  also  suggested  by  different  instability  modes  that  appear  in  the  compressible  flows.  Studies  similar  to  the 
incompressible  studies  shown  in  [63]  appear  to  be  both  interesting  and  very  relevant  to  issues  of  compressible 
turbulence  modeling. 

Experimental  differences  in  numerical  and  laboratory  data.  There  appears  to  be  a  unanimous  agreement 
in  the  DNS  that,  with  increasing  gradient  or  convective  Mach  number,  an  increases  while  022  and  ai2  both 
decrease.  This  trend  does  not  appear  in  the  laboratory  mixing  layer  experiments.  At  low  Me,  {uu)  appears 
to  decrease  more  rapidly,  with  increasing  compressibility,  than  {vv).  This  implies  that  on  decreases  with 
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compressibility  in  contradistinction  to  DNS  results.  It  is  crucial  to  understand  the  cause  for  the  differences 
between  the  DNS  and  the  laboratory  flows  of  the  {uu)  Reynolds  stress  behavior;  turbulence  models  are 
based  on  intuition  gleaned  from  DNS  data  but  are  used  to  predict  engineering  flows  that  are  more  similar 
to  laboratory  flows.  The  present  computational  results  are  closer  to  either  of  the  laboratory  experiments  of 
[53]  or  [54]  than  the  laboratory  experiments  of  [53]  or  [54]  are  to  each  other.  Earnest  speculation  on  the 
source  of  the  differences  in  the  two  experiments  and  facilities  is  required. 


7.  Summary  and  conclusions.  Progress  towards  the  development  of  a  compressible  turbulence  clo¬ 
sure,  starting  at  the  level  of  second-order  moment  equations,  has  been  described.  Modeling  from  the  second- 
order  level  accommodates  important  structural  changes  that  appear  in  the  anisotropy  and  are  a  feature  and 
a  function  of  compressibility.  In  the  second-order  moment  equations  the  compressible  contributions  to  the 
pressxue-strain  covariance  have  been  obtained.  The  pressure-strain  has  been  closed  by  assuming  that,  as 
is  consistent  with  the  weakly  compressible  limit,  it  can  be  modeled  as  a  tensor  polynomial  linear  in  the 
Reynolds  stresses.  The  difference  from  the  incompressible  case  is  that  the  trace  of  the  compressible  strain  is 
not  zero;  it  is  set  equal  to  the  pressure-dilatation  for  which  models  exist.  The  compressible  pressure-strain 
closure  features  a  dependence  on  several  turbulence  descriptors:  the  turbulent  Mach  number,  the  relative 
strain,  the  gradient  Mach  number,  the  production  and  the  dissipation.  As  a  consequence  the  coefficients  in 
the  compressible  pressure-strain  closure  are  not  constants  but  functions  of  the  compressible  parameters. 

In  addition  to  devising  a  closure  for  the  pressure-strain,  a  closure  for  the  mean  acceleration/mass  flux  terms 
appearing  in  the  Reynolds  stress  equations  has  also  been  developed.  For  the  flows  studied  in  this  article, 
which  are  limited  to  the  flows  for  which  relevant  experimental  data  are  available,  the  acceleration/mass  flux 
moments  are  not  important.  The  mass  flux  terms  will  be  important  in  the  combusting  or  hypersonic  flows 
which  motivated  the  thesis  [37],  of  which  this  article  treats  one  aspect. 

Having  closed  the  compressible  Reynolds  stress  equations  standard  tensor  representation  procedures  have 
been  used  to  produce  a  compressible  algebraic  Reynolds  stress  model  useful  for  flows  near  a  structural 
equilibrium.  A  noteworthy  item  of  this  portion  of  the  work  is  that  the  role  that  the  mean  pressure  and 
mean  density  gradients  plays  on  the  Reynolds  stresses  is  immediately  seen.  As  a  consequence,  it  is  seen 
that  the  baroclinic  dyad  can,  in  the  absence  of  mean  velocity  gradients,  contribute  to  the  anisotropy  of  the 
turbulence.  It  is  also  seen  that  the  mean  bulk  dilatation  contributes  to  the  anisotropy  in  a  way  that  is  quite 
different  from  an  irrotational  strain. 

The  mathematical  results  developed  here  have  been  implemented  in  mixing  layer  computations  spanning  a 
wide  range  of  mean  flow  Mach  numbers.  In  this  article  the  discussion  has  been  limited  to  the  compressible 
miying  layer  for  which  a  sizable  amount  of  literature  exists.  The  calculations  presented  have  been  organized 
along  two  themes:  1)  an  investigation  of  the  effects  of  compressibility  as  related  to  the  compressible  pressure- 
strain  and  the  Reynolds  stresses  and  2)  a  validation  of  the  algebraic  Reynolds  stress  model  predictions. 

The  computations  with  the  compressible  pressure-strain  indicate  that  the  present  modeling  [which  does  not 
have  undefined  tunable  constants]  produces  precisely  the  behavior  seen  in  the  DNS  of  [14]  and  [15];  there  is  a 
commensurate  reduction,  with  increasing  convective  Mach  number,  of  all  the  components  of  the  compressible 
pressure-strain  tensor.  The  changes  in  the  pressure-strain  lead,  as  is  established  by  DNS  of  several  different 
flows  [14,  15,  16],  to  changes  in  the  anisotropy  of  the  Reynolds  stresses.  The  changes  predicted  by  the 
modeling  for  the  normal  and  shear  anisotropies  are  very  consistent,  in  trend,  with  DNS  data  and  especially 
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with  the  DNS  of  the  mixing  layer  of  [14]  —  the  flow  conflguration  most  similar  to  the  one  treated  in  this 
article.  Comparisons  between  the  computational  predictions  and  laboratory  resutls  for  the  Reynolds  stress 
and  their  behavior  with  increasing  compressibility  is  poor.  This  is  due  to  unknown  experimental  issues: 
the  agreement  between  even  the  different  laboratory  experiments  is  poorer  than  the  agreement  between  the 
numerical  and  experimental  results. 

All  the  laboratory  data  do  however  agree  on  the  trend  of  the  mixing  layer  growth  rate  with  increasing 
compressibility:  it  decreases.  The  computational  experiments  conducted  indicate  that  sizable  reductions  in 
the  mixing  layer  growth  rate  accompany  changes  in  the  anisotropy  of  the  turbulence  due  to  the  compressible 
aspects  of  the  pressure-strain  covariance. 
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Appendix  A.  The  Poisson  equation  for  the  fluctuating  pressure  is  obtained  by  taking  the  divergence 
of  the  Navier-Stokes  equations.  Both  instantaneous  and  averaged  forms  of  these  equations  are  necessary  in 
the  derivation.  The  subtraction  of  the  mean  equations  from  the  instantaneous  equations  and  simple  term 
manipulations  involving  the  continuity  relation  provide  the  desired  Poisson  equation  [5]: 
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The  solution  for  the  fluctuating  pressure  can  be  determined  by  using  Green’s  functions  or  Fourier  transforms. 
With  respect  to  the  incompressible  counterpart,  the  compressible  Poisson  equation  is  severely  complicated. 
An  equivalent  form  which  resembles  more  to  the  incompressible  equation  can  be  obtained  as: 
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The  last  term  on  the  right  hand  side  (RHS)  converts,  when  the  fluctuating  pressure  solution  is  used  to  deter¬ 
mine  the  pressure-strain  correlation,  into  a  viscous  interaction  tensor  whose  trace  is  equal  to  the  dilatational 
dissipation.  In  the  absence  of  walls  this  term  can  be  omitted  in  high  Reynolds  number  flows.  The  acoustic 
term  (the  first  term  on  the  RHS)  is  difficult  to  be  taken  into  account  in  the  present  analysis.  If  it  is  assumed 
that  the  pressure  ffuctuations  are  caused  by  turbulence  only  then  this  term  is  negligible.  The  condition 
p^/p  1  allows  the  neglect  of  the  second  term  on  the  RHS.  The  remaining  terms  are  the  return-to-isotropy 
and  the  rapid  part  for  which  known  modeling  principles  can  be  applied  in  the  limit  of  low  convective  Mach 
numbers. 


Appendix  B.  It  well  known  that  for  linear  pressure  strain  forms  it  is  impossible  to  satisfy  realizability 
conditions  -  the  requirement  that  the  eigenvalues  of  the  Reynolds  stress  tensor  remain  positive.  Satisfying 
realizability  is  a  very  practical  computationally  stabilizing  requirement.  Our  experience  with  computations 
in  complex  flows  indicates  that  realizability  is  very  useful.  The  model  is  now  made  (weakly)  realizable 
following  methods  suggested  by  Schumann  [64]  and  Lumley  [34]  and  detailed  by  Shih  and  Shabbir  [65].  Let 
F  =  1  +  27J///8  +  9/7/4  is  a  parameter  involving  the  second  invariant  II  =  -^aijaji  and  third  invariant 
III  —  —laijajkaki  of  the  Reynolds  stress  anisotropy  tensor.  Then  the  following  asymptotic  behavior  for  the 
pressure  strain-model  ensures  that  realizability  is  satisfied: 

2 
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o 

(82)  ^^Xpege  ^O  as  F  ^  0 

where  the  index  e  indicates  that  the  relations  are  written  in  the  principal  axes  of  {u^/uj).  The  computa¬ 
tional  form  of  the  model  with  the  additional  parameters  necessary  to  enforce  realizability  so  necessary  for 
computational  stability  is 
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(83)  [1  -  C4  -  2(i2][aipfJpj({'u))  -  fiip((u))apj]  +  ^d2Spp{{u))aij  BrF^’’ 

with  C3  =  (5  —  9(72) /II  and  C4  =  (l  +  7(72)/ll.  The  value  for  the  constant  C2  will  be  the  same  as  in 
the  incompressible  model  to  preserve  consistency  in  the  zero  Mach  number  limit,  that  is  C2  —  0.45.  The 
parameters  are  ar  =  0.1,  /3r  =  0.5,  Ar  =  min(F“"'',  0.1“"’')  and  Br  =  min(jF“^’’,  0.1“^’’). 

Using  this  modified  form  of  the  pressure-strain  model  the  ARSM  coefficients  become:  61  =  |  —  BrF^’-{^  + 
|di),  62  =  bs  =  1  -  BrF^-{l  -C3  +  2d2),  64  =  1  -  BrFl^-{l  -  (74  -  2^2),  and 
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Fig.  1.  Comparisons  between  ARSM  and  SOM  calculations  for  the  mixing  layer,  Me  =  1.07;  Top:  streamwise  normal 
Reynolds  stress;  Bottom:  Reynolds  shear  stress,  fx.  =  1.2/ 
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Fig.  2.  Comparisons  between  ARSM  and  SOM  calculations  for  the  mixing  layer.  Me  =  1.07;  Top:  crosstream  normal 


Reynolds  stress;  Bottom:  turbulent  kinetic  energy  production  over  dissipation,  ix.  =  1.2/ 
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Fig.  3.  Influence  of  the  compressibility  correction  in  SOM  calculations  for  the  mixing  layer,  Me  —  1.07;  Top:  the  voHicity 


thickness;  Bottom:  components  of  the  pressure- strain  tensor.  ^  =  1.2/ 
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Fig.  7.  Effect  of  compressibility  on  normalized  thickness  growth  rate,  k  .  e  calculations  of  a  mixing  layer,  fx  =  2/ 
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Fig,  8.  Effect  of  compressibility  on  normalized  thickness  growth  rate,  ARSM  calculations  of  a  mixing  layer. 
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